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EXECUTIVE  SUMMARY 


1 . Introduction 

This  is  the  final  report  for  Phase  I (first  year)  of  a two  year  effort 
on  Contract  D0T-0S-40093,  "Improved  Wheel  and  Rail  Performance  via  Control 
of  Contact  Stress."  The  general  state  of  the  art  prior  to  the  beginning 
of  the  project  has  been  summarized  in  the  report  "A  Review  of  Rail -Wheel 
Contact  Stress  Problems,"  by  B.  Paul,  FRA-OR&D  76-141,  PB251 238IAS , April  1975 
The  present  report  gives  the  detailed  mathematical  theory  of  a new  approach 
to  the  higherto  unsolved  problem  of  finding  the  stresses  between  two  closely 
fitted  or  "worn-in"  metallic  surfaces,  such  as  a moderately  worn  wheel  and  rail 
Before  applying  the  general  technique  to  the  wheel -rail  problem  it  is  essen- 
tial to  check  its  validity  with  simpler  geometries  such  as  closely  fitting 
cylinders  and  spheres,  where  previous  experimental  and  approximate 
analytical  solutions  exist. 

2.  Problem  Statement 

The  overall  objective  of  the  contact  is  to  generate  a method  for 
calculating  the  contact  stresses  between  arbitrarily  profiled  wheel  and  rails. 
In  this  report  a general  approach  to  the  problem  is  formulated,  and  applied 
to  two  specific  geometries:  (a)  a cylinder  pressed  against  a closely  fitted 

cylindrical  seat  and  (b)  a sphere  pressed  against  a closely  fitted  spherical 
socket.  In  addition,  the  stress  concentrations  induced  by  the  presence 
of  a small  defect  such  as  a corrosion  pit  are  calculated  for  the  case  of 
a sphere. 

3.  Results  Achieved 

Since  worn  wheels  and  rails  contact  conformally,  the  existing  contact 
stress  theories  for  nonconformal  contact  are  not  adequate.  In  this  report  a 
general  numerical  method  of  solution  for  three  dimensional,  frictionless, 
conformal , elastic  contact  problems  is  presented  for  the  first  time.  The 
method  is  used  to  analyze  the  conformal  contact  of  a sphere  indenting  a 
spherical  seat  and  a cylinder  indenting  a spherical  seat.  The  results  of  the 
sphere-spherical  seat  problem  compared  well  with  experimental  data  and  are 
significantly  more  accurate  than  those  of  a previously  published  attempt  to 
solve  the  problem.  Results  of  the  cylinder-cylindrical  seat  problem  were  in 
close  agreement  to  a known  approximate  solution  of  this  problem  and  agree  well 
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with  an  existing  photoelastic  experiment.  For  both  analyses,  results  compared 
favorably  with  Hertzian  theory  for  the  limiting  case  of  small  contact  regions. 

A method  is  given  for  defining  the  boundaries  of  the  large  contact  regions, 
and  for  solving  the  associated  governing  singular  integral  equation  of  the  first 
kind.  A general  iterative  procedure  is  developed  which  converges  to  the 
true  three-dimensional  contact  region. 

In  addition  the  solution  to  a non-Hertzian  contact  problem  with  a multiply 
connected  contact  region  is  solved;  namely,  the  case  of  two  spheres  in  con- 
tact where  one  of  them  has  a surface  defect  or  pit.  The  method  developed  was 
capable  of  detecting  extremely  steep  gradients  in  stress  at  the  defect. 

4.  Utilization  of  Results 

Better  understanding  of  the  contact  stress  distribution  at  the  interface 
of  wheel  and  rail  could  lead  to  substantial  advances  in  the  solution  of  several 
key  problems  in  railroad  technology.  Examples  include  wheel  screech,  flange 
impact  noise,  wheel  and  track  wear  and  fatigue  failures,  deterioration  of 
ride  quality  and  possible  derailment  due  to  lateral  and  longitudinal  slippage, 
increase  of  headway  (and  loss  of  economic  capacity)  due  to  adhesion  limits 
in  braking  and  acceleration. 

The  results  of  the  research  has  potential  for  wheel  and  rail  designers, 
and  those  doing  research  and  development  work  in  the  areas  of  wheel  and 
rail  failures,  rail-car  dynamics,  ride-comfort,  and  safety. 

5.  Concl usions 

The  objectives  set  for  Phase  I of  the  research  have  been  achieved.  In 
addition  to  the  generation  of  a comprehensive  survey  report  on  wheel -rail 
contact  stress  phenomena,  the  work  reported  on  herein  successfully  tested  the 
validity  of  a new  method  for  finding  contact  stresses  between  closely  fitted 
curved  surfaces  such  as  cylinders  and  spheres.  This  work  is  a necessary 
precursor  to  the  solution  of  the  more  complicated  geometrical  configuration 
of  wheel-rail  interfaces,  which  is  the  subject  of  Phase  I of  this  research 
project. 
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1.  INTRODUCTION 


The  very  large  contact  stresses  which  exist  between  rails  and  conven- 
tional wheels  may  be  calculated  by  Hertz's  analysis  when  the  wheels  are  new, 
and  the  area  of  contact  is  small.  However,  when  the  wheels  are  worn,  or  are 
initially  fabricated  with  so-called  preworn  profiles,  the  contact  area  will 
be  too  large  for  the  Hertzian  analysis  to  be  valid.  In  fact,  for  this 
latter  case  of  so-called  conformal  elastic  contact,  there  is  no  currently 
available  method  for  accurately  predicting  contact  stresses.  We  have 
therefore  undertaken  the  task  of  developing  general  methods  for  the  deter- 
mination of  contact  regions,  surface  deformations,  approach^,  and  interfacial 
pressures  in  conformal  (i.e.  closely  fitting)  elastic  bodies.  In  this 
work  we  report  upon  the  numerical  method  developed  to  date,  and  show  how 
it  may  be  applied  to  the  case  of  conformal  cylinders,  or  spheres. 

Contact  problems  can  be  classified  into  the  following  two  categories: 
i)  Problems  where  one  body  is  elastic  and  the  other  is  rigid 
ii)  Problems  involving  two  elastic  bodies 
In  the  first  class  of  problems,  termed  "punch"  problems,  the 


The  approach"  is  defined  in  contact  mechanics  as  the  displacement  of 
a point  in  one  body  relative  to  a point  on  the  other  body,  where  both  points 
are  far  removed  from  the  region  of  contact. 
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The  essentials  of  this  work  constitute  the  Ph.D.  dissertation  of 
W.  Woodward  at  the  University  of  Pennsylvania,  1976. 
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contact  region  is  known  a-priori.  In  the  second  class,  termed  elastic 
contact  problems,  the  contact  region  is  initially  unknown  and  must  be 
determined.  The  first  widely  acclaimed  solution  to  a contact  problem 
was  that  published  by  H.  Hertz  [1881]  involving  the  elastic  contact  of 
frictionless  bodies  with  quadratic  surfaces.  Hertz's  solution  is 
centered  around  the  assumption  that  the  dimensions  of  the  contact 
region  are  small  compared  to  the  radii  of  curvature  of  the  bodies. 
Problems  for  which  this  assumption  is  valid  are  termed  "nonconformal " 
or  "counterformal " contact  problems.  Most  solutions  that  have  been 
found  to  date  are  of  this  type.  In  contrast,  problems  not  restricted 
to  this  assumption  are  termed  "conformal."  Following  Hertz,  solutions 
to  punch  problems  were  analyzed  by  several  Russian  authors  such  as 
Muskhelishvili  [1963],  For  detailed  accounts  of  these  problems  the 
reader  is  referred  to  the  excellent  reviews  of  this  work  by  L.  A. 

Galin  [1961]  and  A.  I.  Lure  [1964].  Recently,  elastic  contact  prob- 
lems involving  friction  and  dynamics  have  also  been  analyzed.  In  a 
recent  publication,  Kalker  [1975]  categorizes  the  solutions  to  date 
and  identifies  the  areas  within  contact  mechanics  which  need 
investigation.  His  comparison  reveals  that  the  areas  involving 
friction,  plasticity,  visco-elasticity  and  large  deformations  are  in 
most  need  of  study.  Kalker  does  not  review  the  analysis  of  conformal 
contact  problems  in  his  survey.  It  should  be  noted  that  the  uniqueness 
theorem  of  Kirchoff  (not  intended  for  contact  problems)  was  extended 
only  recently  to  include  general  frictionless,  elastic  contact  problems 
by  Kalker  [1971] . 
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This  dissertation  centers  on  providing  a general  numerical 
method  of  solution  to  conformal  frictionless  contact  problems.  The 
literature  in  the  areas  mentioned  above  is  far  too  comprehensive  to 
review  in  this  brief  introduction,  instead  the  interested  reader  is 
referred  to  the  aforementioned  surveys.  The  remainder  of  this  dis- 
cussion will  be  devoted  to  a more  detailed  discussion  of  the  existing 
numerical  solutions  and  the  limited  literature  on  conformal  contact 
probl ems . 

With  the  advent  of  the  digital  computer  several  numerical 
techniques  have  been  developed  to  analyze  a more  general  class  of  con- 
tact problems.  Conry  and  Seireg  [1971]  have  examined  elastic  contact 
in  terms  of  a linear  programming  model.  Their  method  is  general  in 
scope,  however,  the  only  examples  which  were  analyzed  were  Hertzian  or 
one  dimensional  beam  problems. 

Kalker  and  van  Randen  [1972]  derived  a variational  principle 
for  both  linear  and  non-linear  elastic  contact  problems.  For  the  case 
of  linear  elasticity  the  principle  takes  the  form  of  an  infinite 
dimensional  convex  quadratic  programming  problem.  They  successfully 
solved  both  a Hertzian  and  non-Hertzian  problem.  It  was  concluded  that 
the  solution  yielded  accurate  values  of  approach,  maximum  pressure  and 
applied  force;  however,  the  actual  contact  area  was  not  determined  with 
great  accuracy. 

Finite  element  techniques  have  also  been  adapted  to  solve 
contact  problems  by  Chan  and  Tuba  [1971]  and  more  recently  by  Chaud, 
Haug  and  Rim  [1974].  Both  methods  are  general  in  that  they  are 
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reported  to  be  able  to  handle  problems  which  fall  into  the  domain  of 
finite  element  analysis  such  as  analyzing  non-isotropic, 
non-homogeneous  media  or  problems  with  plasticity  and  creep,  however, 
both  works  report  only  examples  which  are  composed  of  isotropic 
materials  stressed  within  the  range  of  linear  elasticity. 

Tuba  and  Chan  compare  their  computed  results  to  photo  elastic 
studies  and  concluded  that  trends  were  identical  but  the  results  lacked 
close  agreement.  Chaud  et  al,  analyzed  the  non-Hertzian  problem  of  a 
human  knee  joint  and  the  contact  between  two  half  spaces  where  one 
half  space  has  three  bumps  on  the  surface.  The  contact  area  in  the 
latter  case  found  in  photo  elastic  studies  had  good  general  agreement 
with  their  computed  results. 

A general  method  of  solution  of  non-Hertzian,  non-conformal 
elastic  contact  problems  was  developed  by  Singh  and  Paul  [1974].  They 
considered  the  classical  contact  criterion  (which  includes  that  of 
Hertz)  for  arbitrary  surface  geometries.  In  order  to  solve  the 
governing  integral  equation  of  the  first  kind,  which  belongs  to  the 
class  of  "ill  posed"  or  "Hadamard  incorrect"  problems,  they  introduced 
three  different  numerical  schemes.  The  first  "simply-discretized 
method"  was  found  to  be  relatively  unstable  for  the  particular  problems 
they  investigated.  In  order  to  overcome  this  difficulty,  Singh  and 
Paul  [1973-74]  introduced  two  other  methods  of  solving  ill  posed 
integral  equations,  called  the  "Redundant  Field  Point  method"  and 
the  "Functional  Regularization  method";  the  latter  of  which  is  based 
on  Tychonov's  regularization  procedure. 
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In  contrast  to  counterformal  problems  are  those  of  the 
"conformal"  type;  i.e.,  those  where  the  dimensions  of  the  contact 
region  can  be  large  compared  to  the  smallest  radii  of  curvature  of  the 
contacting  bodies.  Relatively  few  solutions  to  conformal  contact 
problems  have  been  published.  A brief  summary  of  all  elastic  conformal 
problems,  known  to  the  author,  follows. 

An  elastic  sphere  indenting  an  elastic  seat  has  been  solved  by 
Goodman  and  Keer  [19651  . They  present  results  for  the  half  angles  of 
contact  up  to  20  degrees  and  provide  experimental  results  which 
generally  agree  with  their  solutions.  Improvements  to  the  Hertzian 
theory  are  discussed — in  particular,  the  problem  which  arises  when  one 
tries  to  include  terms  of  higher  order  than  those  used  in  the  "half 
space  solution"  (of  Boussinesq)  which  is  fundamental  to  the  Hertzian 
solution.  It  is  noted  that  there  are  higher  order  terms  in  the  exact 
formulation  of  the  sphere  problem  which  do  not  appear  in  the  formula- 
tion if  the  half  space  assumption  is  used  without  truncating  terms. 
These  terms  are  particular  to  the  spherical  geometry.  Goodman  and  Keer 
justify  their  extension  of  the  Hertzian  theory  through  analysis  of 
these  second  order  terms. 

The  conformal  contact  of  an  elastic  cylinder  indenting  a 
cylindrical  seat  was  first  analyzed  by  Sjtaerman  [1940]  and  more 
recently  by  Persson  [1964].  Sjtaerman  and  Persson  derived  the  iden- 
tical "contact  criterion"  but  both  proceeded  in  different  ways  to  solve 
the  equation.  Sjtaerman  formulated  the  displacements  in  terms  of 
integrals  of  the  influence  functions  and  used  finite  difference 
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techniques  to  solve  the  resulting  integral  equation  for  the  unknown 
pressure  distribution.  On  the  other  hand,  Persson  assumed  the  contact 
region  to  be  cylindrical  in  shape  and  formulated  the  criterion  as  an 
integro-differential  equation  from  which  he  found  analytic  expressions 
for  the  pressure  field.  The  earlier  solution  of  Sjtaerman  appears  to 
be  less  accurate,  possibly  because  he  published  before  the  digital 
computer  was  invented  and  may  have  been  forced  to  use  a too  crude 
finite  difference  mesh. 

Recently,  a number  of  problems  involving  a disc  in  an  infinite 
plate  under  tension  have  been  solved  by  finite  element  techniques. 

Chan  and  Tuba  [1971]  analyzed  a plate  under  tension  with  a shrink  fit 
disc  located  in  the  center.  They  present  results  which  show  good 
agreement  between  their  computed  values  of  ci rcumferential  stress  and 
the  exact  solution,  however,  there  is  a larger  discrepancy  between  the 
computed  value  of  compressive  stress  and  the  exact  solution.  In  fact 
the  compressive  stress  on  each  body  for  any  one  angle  does  not  in 
general  agree. 

Chaud  et  al  [1974]  have  analyzed  the  problem  of  a plate  under 
tension  with  either  a loose  or  full  inclusion.  They  show  good  agree- 
ment between  their  predicted  contact  stress  and  experimental  results 
for  a contact  angle  of  20  degrees. 

The  goal  of  this  research  is  to  develop  a general  method  of 
analysis  for  frictionless,  conformal  contact  problems.  In  particular, 
the  method  developed  is  to  be  used  in  future  research  on  the  analysis 
of  interfacial  contact  stresses  between  a railway  wheel  and  rail.  The 
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concepts  of  Hertz's  classical  geometric  formulation  of  the  contact 
criterion  are  extended  to  non-planar  surfaces  resulting  in  a singular 
integral  equation  of  the  first  kind.  The  solution  is  dependent  on 
identifying  the  influence  functions  for  surface  point  displacements 
in  the  region  of  the  contact  area.  A numerical  appraoch  to  generating 
influence  functions  is  developed  and  the  accuracy  of  the  generated 
functions  is  shown  to  be  good  for  those  cases  where  analytic  solutions 
are  known.  The  simply  discretized  method  of  Singh  and  Paul  [1973]  was 
used  for  solution  of  the  integral  equation.  The  solutions  using  this 
formulation  were  compared  to  Hertz's  solutions  for  limiting  cases  in- 
volving small  contact  regions. 

The  results  of  the  present  general  method  are  compared  to 
available  analytic  solutions  to  specific  problems  involving  the  contact 
of  an  elastic  cylinder  in  seat  and  a sphere  in  seat. 

In  addition,  a solution  was  found  to  a noncomformal  problem 
with  a multiply  connected  contact  region.  The  proper  boundary 
iteration  which  is  necessary  to  arrive  at  a unique  solution  is 
developed  and  discussed.  The  specific  example  analyzed  is  a pitted 
sphere  in  contact  with  a sphere.  The  significance  of  pit  geometries 
on  the  contact  stress  are  illustrated. 

In  summary,  the  main  contributions  of  this  disseration  are: 

1.  A general  numerical  method  for  solution  of  frictionless  conformal 
elastic  contact  problems  is  presented 
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2.  Numerical  influence  functions  needed  for  solutions  of  problems 
with  cylindrical  and  spherical  surface  geometries  were  generated 
and  their  accuracy  was  verified  by  comparison  to  exact  analytic 
solutions  when  they  existed 

3.  In  its  application  to  the  specific  problems  of  a sphere  indenting 
a spherical  seat  and  a cylinder  indenting  a cylindrical  seat,  it 
was  shown  that  this  method  produces  accurate  values  of  contact 
pressure  approach,  displacements,  strains  and  applied  force 

4.  The  problem  of  a pitted  sphere  indenting  a sphere  was  solved  for 
the  first  time  and  the  appropriate  boundary  iteration  for 
multiply  connected  contact  regions  was  established 

Chapter  2 contains  a brief  review  of  the  previous  non- 
conformal  methods  of  solution  presented  by  Singh  and  Paul  [1973-74]. 
The  conformal  contact  theory  which  is  the  basis  of  this  work  is 
formulated  in  chapter  3.  The  generation  of  influence  functions  which 
are  necessary  to  the  solution  of  the  sphere  and  cylinder  problems  is 
discussed  in  chapter  4.  Chapters  5 and  6 contain  solutions  to  the 
examples  of  a sphere  indenting  a conformal  spherical  seat  and  a cyl- 
inder indenting  a conformal  cylindrical  seat  respectively.  A contact 
problem  involving  a multiply  connected  contact  region  is  solved  in 
chapter  7.  The  conclusions  of  this  work  are  presented  in  chapter  8. 


2.  FORMULATION  AND  SOLUTIONS  FOR  NON-CONFORMAL  CONTACT  PROBLEMS 


2.1  Introduction 

The  basic  equation  for  non-conformal  contact  theory  is 
developed  in  this  chapter  along  with  a method  of  solution.  Hertz 
[18811  has  shown  that  the  governing  equation  is  an  integral  equation 
of  the  first  kind.  Hertz  found  an  analytical  solution  to  this 
equation  for  the  special  case  where  the  surfaces  may  be  modelled  as 
locally  quadratic;  however,  the  integral  equation  itself  applies  to 
any  non-conformal  contact  problem  and  has  been  solved  by  Singh  and 
Paul  [1974]  for  non-conformal,  non-Hertzian  contact  problems.  The 
method  of  solution  outlined  here  is  that  developed  by  Singh  and  Paul 
and  is  termed  the  "Simply-Discretized"  or  "S.D."  method.  They  proved 
that  the  S.D.  method  can  become  unstable;  they  applied  a stabilizing 
technique  termed  the  "Functional  Regularization"  or  "F.R."  method, 
when  the  S.D.  method  proved  unstable  and  successfully  solved  several 
problems.  The  "Functional  Regularization"  method  is  also  summarized 
in  this  chapter. 

In  addition  to  presenting  the  basic  integral  equation  for 
non-conformal  contact,  the  material  in  this  chapter  introduces  the 
concepts  of  contact  theory  which  will  be  used  in  the  development  of 
the  governing  integral  equation  for  conformal  contact  in  chapter  3. 
Furthermore,  the  "Simply-Discretized"  method  of  solution  will  be  used 
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in  solving  the  conformal  contact  problems  presented  in  chapters  5 and 
6 and  in  the  non-conformal  contact  problem  with  a multiply  connected 
contact  region  in  chapter  7. 

2.2  The  Governing  Equation  for  Non-conformal  Contact  Theory 

Consider  two  frictionless  non-conforming  bodies  initially  in 
contact  at  a single  point.  Loading  each  body  such  that  the  resultant 
force  acts  through  the  initial  point  of  contact  produces  deformations 
in  the  neighborhood  of  the  initial  point  of  contact.  The  area  of 
contact  between  the  bodies  will  increase  from  a single  point  to  a 
finite  area.  In  non-conformal  contact  theory  it  is  assumed  that  the 
dimensions  of  the  contact  area  are  small  compared  to  the  local  radii 
of  curvature  of  the  two  contacting  surfaces.  After  deformation,  the 
two  bodies  undergo  a localized  elastic  deformation  and  a rigid  body 
displacement.  The  rigid  body  displacement  is  referred  to  as  the 
"approach"  in  contact  mechanics. 

In  general  the  geometry  of  the  surfaces  before  deformation 
and  the  applied  thrust  is  known,  while  the  actual  contact  area,  the 
pressure  distribution  within  this  area,  and  the  approach  6 are  unknown. 
The  governing  equation  for  non-conformal  contact  relates  the  approach, 
the  contact  area,  the  surface  geometry  and  the  interfacial  pressure 
distribution. 

Consider  body  1 and  body  2 initially  in  non-conformal  contact 
at  a point  0.  (Fig.  2.1)  Let  a right-hand  cartesian  coordinate  sys- 

AAA 

tern  (x.j,  y^,  z.j)  be  constructed  such  that  the  x-j  - y-|  plane  lies 

A 

tangent  to  body  1 at  point  0 and  has  point  0 as  its  origin.  Let  z-j  be 
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Fig.  2.1.  Coordinate  systems  for  non-conformal  contact 
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the  unit  inward  normal  to  body  1.  Also  define  a left-handed 

AAA  A A A A 

coordinate  system  (x2,  y2,  z2)  such  that  x2  = x] » ^2  ~ yl ’ and 

A A 

z2  = - z 1 . Now  examine  a cross  section  of  the  contacting  bodies  through 
the  point  0 as  shown  in  figure  2.2.  It  is  assumed  that  points  on  the 
surfaces  of  bodies  1 and  2 which  merge  after  deformation  are  located 
at  the  same  (x,  y)  coordinates.  The  displacements  due  to  the  contact 
phenomena  of  a point  A on  body  1 and  B on  body  2 which  merge  after 
deformation,  will  be  examined  closely  in  the  following  paragraphs. 

Consider  the  change  in  the  position  of  point  A on  body  1 
after  a load  F is  applied.  Due  to  the  elastic  deformation  of  the 
surface,  point  A moves  the  distance  A 'A"  in  the  z-j  direction.  This 
elastic  deformation  will  be  labeled  w-j . Also  point  A moves  a distance 
AA'  due  to  rigid  body  motion,  labeled  A Similarly  on  body  2 a point 
B displaces  an  amount  B'B"  or  w2  due  to  elastic  deformation  and 
displaces  from  B to  B1  due  to  a rigid  body  motion  A Therefore,  con- 
sidering the  total  motion  of  points  A and  B,  point  A moves  a total 
distance  w-j  + and  point  B moves  w 2 + A^. 

Having  traced  the  motion  of  points  on  the  surfaces  of  the 
contacting  bodies,  these  motions  may  be  related  to  the  surface 
geometry  to  obtain  a necessary  condition  for  contact  of  the  two  sur- 
faces. Note  that  points  A and  B are  initially  separated  by  a distance 
equal  to  f-j  (x,  y)  + f2  (x,  y)  where  f-j  (x,  y)  and  f2  (x,  y)  are 
termed  the  "profile  functions"  of  the  two  surfaces,  f^  (x,  y) 
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Fig.  2.2.  Kinematics  of  surface  point  displacements  in  non- 
conformal  contact. 
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represents  the  coordinate  of  a point  (x,  y,  z-j ) on  the  surface  of 
body  1,  while  f^  (x,  y)  is  the  z^  coordinate  of  a point  (x,  y,  z^)  on 
the  surface  of  body  2.  After  deformation  the  initial  separation 
changes.  Consider  the  final  separation,  S,  of  points  A and  B after 
deformation 

s*  + (w,+a>)*(Wi+Ai)  (2.i) 


Define  a function  f (x,  y)  and  a scalar  6 such  that 

f (*,$)-  +4  (x-v)  <2-2) 

and 

<J=  (2.3) 

The  scalar  6 is  termed  the  "approach"  and  physically  represents  the 
distance  that  points  on  one  body  move  parallel  to  the  z-^axis  towards 
points  on  the  other  body  due  to  rigid  body  movement.  Rewriting 
equation  (2-1)  by  substituting  equations  (2-2)  and  (2-3), 

S-  ) + (w,  t wz ) - S (2.4) 

By  assumption  in  non-conformal  contact  theory,  the  contact  area  is 
small  compared  to  the  local  radii  of  curvature  of  the  contacting 
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bodies.  Therefore,  it  is  appropriate  to  replace  these  displacements 
with  the  displacement  field  of  an  elastic  half  space  due  to  some 


where  and  v.  are  the  Young's  modulus  and  Poisson's  ratio  for 


body  i.  It  is  reasonable  to  neglect  the  displacements  in  the  (xj  - yp 

/s  /\ 

plane  and  assume  that  points  A and  B,  located  at  equal  x^»y^  coordi- 
nates, will  merge  after  deformation  since  the  displacements  in  the 
(x-j-  y|)  plane  on  each  body  are  nearly  equal  in  magnitude  and  direction. 


The  pressure  distribution  in  equation  (2.5)  is  over  the  contact 
region  ft..  on  body  i,  however  ft-j  = 9,^.  Furthermore,  since  the  radii  of 
curvature  of  the  bodies  are  large  compared  to  the  dimensions  of  the 
contact  area  and  the  radius  of  curvature  of  the  contact  patch  is  still 
larger  owing  to  the  non-conformal  nature  of  the  contact,  the  contact 
area  may  be  represented  by  ft,  the  projection  of  ft-j  onto  the  x^  - y^ 
plane.  Substituting  equation  (2.5)  into  equation  (2.4)  the  separation 


pressure  distribution  p (x,  y)  over  an  area  ft..  Hence, 


2 


(2.5) 


t 
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becomes: 


3 - f(x,  3)  t 


k f 

J J fr-x'JViM')2 

ob&  4m 


(2.6) 


where 


n E,  n 


If  U is  known  then  equation  (2.6)  is  an  integral  equation  of 
the  first  kind  and  is  the  governing  equation  for  non-conformal  contact. 
The  separation  of  two  points  within  the  contact  region  has  to  be  zero 
while  the  pressure  has  to  be  positive,  i.e.,  the  bodies  can  only  exert 
compressive  forces  on  one  another  within  the  contact  region.  Outside 
the  contact  region,  S must  be  positive  while  the  pressure  distribution 
must  equal  zero. 

These  boundary  conditions  may  be  summarized  as  follows: 


S - 0 IS/SIDE  J1 

J >0  outside  n 

<P( x,  i) ! >0  iusipe  sx 

<P(X,  1)  - O 


(2.7a) 

(2.7b) 


(2.7c) 


OUTSIDE 


(2.7d) 
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The  solution  of  a non-conformal  contact  problem  for  a given  thrust 
requires  the  determination  of  the  unknown  contact  region  n,  the 
pressure  distribution  p (x,  y),  and  the  approach  6.  These  quantities 
must  all  satisfy  equation  (2.6)  and  the  boundary  conditions  (2.7). 
General  numerical  solutions  to  equation  (2-6)  are  explained  in 
sections  2.3  and  2.4. 

2.3  The  "Simply-Discretized"  Method  of  Solution  of  Singh  and  Paul 

The  general  non-conformal  contact  problem  as  posed  in  equation 
(2-6)  and  boundary  conditions  (2.7)  has  been  solved  by  Singh  and  Paul 
by  the  "Simply-Discretized"  method.  This  method  is  a semi  inverse 
solution  which  will  be  explained  in  the  following  paragraphs. 

Given  that  two  bodies  of  known  shape  are  brought  into  contact 
with  one  another  and  held  there  by  a force  F,  the  task  remains  to 
locate  the  boundary  of  the  contact  region  fi  and  to  find  the  interfacial 
pressure  distribution  p (x,  y)  and  the  approach  6.  The  "Simply- 
Discretized"  method  is  called  a semi  inverse  method  because  the  contact 
area  is  assumed  to  be  some  logical  "candidate"  region  whereupon  the 
pressure  distribution  and  approach  are  then  found  via  equation  (2.6). 
The  force  F is  then  calculated  from  the  integral  of  the  pressure  dis- 
tribution over  the  area  n.  The  initial  guess  of  the  logical 
"candidate"  contact  region  ft  is  then  checked  to  see  if  the  values  of 
the  separation  satisfy  equation  (2. 7a)  and  (2.7b).  If  they  do  not, 
the  "candidate"  region  can  be  modified  to  better  approximate  the  true 
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contact  region  and  the  solution  procedure  is  repeated.  Now  consider 
the  details  of  this  solution. 


Discretized"  method  is  to  find  an  approximation  to  the  contact 
region  ft.  Singh  and  Paul  used  the  plane  projection  of  an 

"interpenetration  curve"  formed  by  the  intersection  of  the  two 
surfaces  when  one  body  was  allowed  to  mathematically  interpenetrate 
the  other.  Increasing  the  depth  of  interpenetration  would  increase 
the  contact  area.  Each  interpenetration  depth  would  provide  one 
"candidate"  contact  area  corresponding  to  one  loading  F on  the  bodies. 
This  concept  of  interpenetration  is  physically  meaningless  in  a con- 
tact problem  since  the  bodies  can  not  actually  interpenetrate  one 
another;  it  is  only  a method  which  enables  one  to  find  an  approxima- 
tion to  the  contact  area. 


Having  found  a "candidate"  contact  region  the  next  task  is  to 


find  the  pressure  distribution  and  approach  in  equation  (2.6).  A 
"Simply-Discretized"  solution  is  obtained  by  assuming  a piecewise 
constant  pressure  distribution  over  the  area  ft.  Dividing  up  the 
contact  area  ft  into  N cells  and  assuming  the  pressure  p (x,  y)  to  be 
constant  within  each  cell  equation  (2.6)  becomes 


The  first  step  in  solving  equation  (2.6)  by  the  "Simply- 


(2-8) 
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where 


O-v,*) 


TTE 


IT  Ez 


(2.9) 


and 


P.  is  the  constant  pressure  in  cell  i. 


In  equation  (2.8)  there  are  N unkown  P^'s  and  1 unknown  6,  thus  a 
total  of  N + 1 unknowns.  However,  equation  (2.8)  applies  to  every 
point  (x,  y)  inside  n,  thus  it  can  be  written  for  N + 1 "field  points" 
inside  the  contact  region.  Singh  and  Paul  chose  the  centroids  of  the 
N cells  as  N of  the  field  points  and  picked  a last  field  point  at  the 
intersection  of  several  cells.  The  integrals  in  equation  (2.8)  can 
be  evaluated  numerically  to  provide  the  coefficients  for  this  set  of 
linear  equations.  Thus  equation  (2.8)  can  be  expanded  to  the  form 


(2.10) 


where 


P.  is  the  pressure  in  cell  j 

J 

f.j  is  the  initial  separation  of  the  centroid  of  cell  i 
ft.  is  the  area  of  cell  j 

J 

(x. , y.)  are  the  coordinates  of  the  centroid  cell  i 


and 
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Vi  PC 


(2.11) 


where 


d%'  a 


la tj(*e-x')x*  0o-3')* 


(x  , u ) are  the  coordinates  of  the  N + 1 field  point 
oo 

f^  + i is  the  initial  separation  of  the  N + 1 field  point 


Combining  equations  (2.10)  and  (2.11)  into  one  set  of  linear  equations 
in  P..  and  6,  the  unknown  P^  and  6 can  be  found,  in  principle,  by  using 
standard  Gaussian  elimination.  This  is  done  by  first  substracting 
equation  (2.11)  from  equation  (2.10),  thus  eliminating  6 from  all 
equations.  The  new  set  of  N equations  formed  can  be  written  as 


(2.12) 


where 
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Equation  (2.12)  is  a set  of  N linear  equations  in  N unknown  P.j's. 

After  solving  this  set  of  equations  for  the  pressure  field,  the 
approach  6 may  be  found  by  back  substitution  of  these  pressures  into 
equation  (2.11). 

Having  obtained  the  pressure  distribution  and  approach  the 
original  "candidate"  contact  region  ft  can  be  verified  via  the  boundary 
conditions  (2.7a,  b).  The  "candidate"  contact  region  can  be  reduced 
if  the  pressure  is  less  than  zero  at  the  boundary  of  ft  or  it  can  be 
extended  if  the  separation  is  negative  outside  the  boundary  of  ft.  The 
problem  can  be  solved  again  if  necessary  to  find  the  pressure  distri- 
bution in  the  connected  region  ft.  When  the  true  contact  area  is  found 
the  total  force  applied  to  the  bodies  can  be  computed  by  integrating 
the  pressure  over  the  contact  area,  i.e.. 


where 


f =J(P(X,  Dahl  fill 

Jl 


(2.13) 


A.,  is  the-area  of  cell  i. 

In  applying  the  S.D.  method  to  a variety  of  problems,  Singh 
and  Paul  found  that  it  was  numerically  unstable  in  the  general  case. 

For  small  cell  densities  the  solutions  obtained  were  good;  however,  as 
the  cell  density  was  increased,  the  solutions  broke  down.  The  pressure 
distribution  became  very  erratic,  changing  from  positive  values  to 
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negative  values  from  cell  to  cell.  The  details  of  these  results  were 
recorded  by  Singh  [1972]  and  Singh  and  Paul  [1973]  [1974].  One  method 
for  eliminating  the  problem  of  ill -conditioning  is  known  as  the 
"Functional  Regularization"  method  which  will  be  explained  in  section 

2.4. 

2.4  "Functional  Regularization" 

The  "Functional  Regularization"  method  is  a technique  of 
stabilizing  an  ill-conditioned  set  of  linear  equations.  This  method 
was  developed  by  Singh  and  Paul  [1973]  in  order  to  extract  a sensible 
solution  from  the  unstable  "Simply-Discretized"  method  of  solving  non- 
conformal  contact  problems. 

The  "Simply-Discretized"  method  yields  solutions  which  have 
wide  variations  in  the  pressure  fields.  Furthermore,  it  was  observed 
that  small  perturbations  in  the  coefficient  matrix  B of  equation  (2.12) 
produce  completely  different  pressure  fields  which  also  vary  radically 
from  cell  to  cell.  Although  the  determinant  of  the  coefficient  matrix 
was  not  singular,  it  appeared  that  there  were  many  solutions  to  the 
set  of  equations  generated  by  the  "Simply-Discretized"  method  of 
solution.  In  order  to  find  the  correct  solution  vector  of  a large 
number  of  vectors  that  satisfy  the  linear  equation  set,  the  F.R.  method 
exploits  the  "smoothness"  property  of  the  pressure  distribution.  The 
F.R.  method  seeks  to  approximately  solve  the  original  set  of  equations 
while  it  simultaneously  minimizes  the  difference  between  the  pressure 
in  neighboring  cells.  This  is  accomplished  by  introducing  a function  $ 


23 


such  that 

# --  f (Pi-Pt„)x  <2-14> 


Solving  the  equation  set  (2.12)  is  equivalent  to  finding  a solution 

2 

to  equation  (2.15)  which  produces  an  e of  exactly  zero. 


(2.15) 


"Functional  Regularization"  seeks  an  approximate  solution  to  equation 

2 

(2.15)  such  that  both  e and  $ are  small  quantities.  Mathematically, 
this  is  equivalent  to  minimizing  a functional 


Y(Pi)  = f Bcj  Pj  -fc  + V §(Pi)  (216) 

where 

o is  a small  parameter  which  controls  the  influence  of  the 
constraint  function  $.  Minimizing  & implies  finding  the 
vector  P..  such  that  - 0 or  solving  the  equation  set 

pp  = Bit  FC  (2.i7) 

cf  r( 


Bounds  for  the  parameter  u have  been  given  by  Singh  and  Paul  [1973  ]. 


24 

The  solution  to  equation  (2.17)  is  stable  in  that  it  yields 
a smoothly  varying  positive  pressure  distribution  which  is  physically 
realistic.  Singh  and  Paul  [1973-74]  have  shown  that  it  closely  agrees 
with  the  exact  solution  for  several  Hertzian  contact  problems. 


3.  FORMULATION  OF  CONFORMAL  CONTACT  PROBLEM 


3.1  Introduction 

The  development  of  nonconformal  contact  theory  was  based  on 
the  assumption  that  the  dimensions  of  the  contact  area  are  small 
compared  to  the  local  radii  of  curvature  of  the  contacting  surfaces. 
This  assumption  is  no  longer  valid  in  conformal  contact.  By  its 
very  nature,  conformal  contact  can  produce  contact  areas  with 
dimensions  as  large  as  the  radii  of  curvature  of  the  surfaces. 

Because  of  this  assumption  in  the  nonconformal  theory  it  was 
appropriate  to  approximate  the  contacting  surfaces  by  two  elastic 
half  spaces  and  to  use  the  Boussinesq  displacement  function  for  a 
point  load  on  a plane  as  the  influence  function^  necessary  for  the 
calculation  of  the  displacement  field.  In  conformal  contact  theory 
the  contact  region  cannot  be  approximated  by  a plane,  and  alternative 
influence  functions  for  the  surfaces  must  be  found. 

Furthermore,  no  longer  can  the  displacements  tangent  to  the 
surface  be  considered  small  as  was  the  case  in  nonconformal  theory. 

A solution  procedure  must  incorporate  both  the  normal  and  tangential 
displacements  of  the  contacting  surfaces  in  the  solution. 

Presented  in  the  following  sections  of  this  chapter  is  a 

Hhe  "influence  function"  relates  the  surface  tractions  to 
the  displacement  field.  It  may  be  sometimes  referred  to  as  a "Green's 
function. " 
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formulation  for  conformal  contact  theory  which  not  only  incorporates 
the  true  influence  function  for  conformal  surfaces  but  also  accounts 
for  the  surface  displacements  both  normal  and  tangent  to  the  surface. 

3.2  Assumptions  in  Conformal  Contact  Theory 

Consider  two  contacting  bodies,  labeled  1 and  2,  which  have 
closely  conforming  surfaces,  i.e.,  they  exhibit  conformal  contact. 

The  initial  point  of  contact  will  be  labeled  0.  Figure  3.1  represents 
a cross  section  through  0 of  the  two  surfaces.  A coordinate  system 
is  constructed  such  that  z is  the  inward  unit  normal  to  body  1 at  0 

✓N 

and  unit  normal  y lies  in  the  plane  of  the  cross  section  at  90° 
clockwise  of  z.  The  intersection  of  the  surfaces  of  body  1 and  2 with 
planes  through  the  z axis  will  be  termed  the  "contour  curves"  of 
the  respective  surfaces.  The  following  assumptions  will  be  made: 

1.  The  surfaces  are  assumed  to  be  frictionless 

2.  The  line  of  the  applied  load  on  the  bodies  in  contact  passes 
through  0 (fig.  3.1 ) 

Considering  only  frictionless  surfaces  reduces  the  complexity 
of  the  contact  problem  significantly,  yet  it  does  not  destroy  the 
usefulness  of  the  solution.  It  is  often  desirable  to  have  friction- 
less surfaces  in  contact  applications.  For  example,  in  the  situa- 
tions involving  bearing  surfaces,  such  as  ball  bearings  or  ball 
joints,  the  surfaces  are  machined  and  lubricated  to  minimize  surface 
friction.  This  assumption  dictates  that  no  shear  tractions  can  be 
applied  to  the  surface  of  either  body  within  the  contact  region. 
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Fig.  3.1.  Cross  section  of  bodies  in  conformal  contact 


Hence,  the  normal  interfacial  pressures  will  be  the  only  surface 
tractions  allowed  within  the  contact  area. 
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Assumption  (2)  requires  that  the  resultant  applied  force 
passes  through  0.  Mathematically  this  assumption  requires  that 
the  following  relationship  hold: 

S r *(-  fn-)  dA  - o (3.D 

A 

Where  ft..  is  the  contact  area  on  body  i,  r is  a vector  extending  from 

/\ 

the  origin,  0,  to  a point  within  the  contact  area,  n..  is  the  unit 
normal  to  body  i of  the  point  defined  by  r,  and  p is  the  interfacial 
pressure  at  the  point  located  by  r.  This  assumption  is  not  required 
in  the  analysis,  rather  it  is  made  to  simplify  the  analysis.  It 
should  be  possible  to  extend  the  present  analysis  to  include  moments 
and  rigid  body  notations. 

3 . 3 Formulation  of  Conformal  Contact  Criterion 

Consider  two  conforming  frictionless  bodies  in  contact.  Body 
1 will  be  called  the  "indentor"  while  body  2 will  be  termed  the 
"seat."  In  the  undeformed  state  these  bodies  contact  at  a point  0. 
Figure  3.2  represents  a cross  section  through  0 of  typical  conformal 

AAA 

contact  surfaces.  A global  coordinate  system  (x,y,z)  is  constructed 

such  that  the  x-y  plane  is  tangent  to  point  0 on  body  1 with  0 as 

/\ 

its  origin  and  z is  directed  inside  body  1. 

One  of  the  difficulties  in  this  class  of  problems  is  the 


t A ^ 

Ah 

per/N/r/oN  or  outward  normal*, 
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Fig.  3.2  Definition  of  coordinate  systems  in  conformal 
contact  theory. 
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identification  of  points  on  each  body  which  come  into  contact  in  the 
deformed  state.  In  order  to  find  the  sets  of  mating  points  we 
introduce  an  iterative  scheme  termed  the  "point-mating  procedure." 

The  details  of  this  scheme  will  be  discussed  later. 

Consider  two  points,  A on  body  1 and  B on  body  2,  which  come 
into  contact  after  deformation.  It  will  be  initially  assumed  that 
points  A and  B lie  at  equal  distances  along  their  respective  contour 
curves  from  the  initial  point  of  contact.  Figure  3.1  illustrates 
this  concept.  The  distance  between  0 and  A measured  along  the 
contour  curve  of  body  1 is  s^  while  the  distance  between  0 and  B 
measured  along  the  contour  curve  of  body  2 is  s^.  In  order  to  form- 
ulate the  contact  criterion  it  will  be  initially  assumed  that 


(3.2) 


/s 


Let  n-j  and  n^  define  the  outward  unit  normals  to  the  surfaces 
of  bodies  1 and  2 respectively.  Shown  in  figure  3.2  are  the 
normal  vectors  at  points  A and  B.  n^  is  directed  ip  degrees  clockwise 
of  the  z direction  while  n^  is  directed  <f>  degrees  counter  clockwise 


of  the  z direction.  For  extremely  conforming  bodies  n^  ~ -n-j . A 
local  coordinate  system  will  be  assigned  to  each  set  of  points  which 
contact  after  deformation.  This  local  coordinate  system  will  have 
a unit  vector  r defined  as  the  mean  of  n^  and  -n-j , i.e.,  the  angle 


between  r and  z is  defined  by  a measured  counter  clockwise,  where 
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(3.3a) 


A unit  vector  t will  be  defined  as  being  90°  clockwise  of  r in  the 
y-z  plane,  and  a ui  unit  vector  will  be  defined  by 


The  displacement  of  point  A (and  B ) due  to  conformal  contact 
as  illustrated  in  figure  3.3  will  be  traced.  Point  A (B)  undergoes 
a rigid  body  translation  A-j  ( ) which  carries  it  in  the  direction 


of  z (-z)  to  point  A'  (B1).  As  shown  in  figure  3.3,  point  A displaces 


an  amountw^in  the  -n^  direction  (from  A'  to  A")  due  to  elastic 
deformation.  Not  shown  in  figure  3.3  are  the  elastic  displacements 
of  points  A and  B in  the  tangent  plane.  These  displacements 
will  be  initially  neglected  in  the  formulation  of  the  contact 
cri  terion. 


Similarly  point  B displaces  to  B"  due  to  an  elastic  displace- 


ment \^2  in  the  - n ^ direction.  The  original  vector  separation  of 


points  A and  B is  labeled  f and  is  a function  of  the  geometry  of  the 


contour  curves.  The  vector  separation  S,  between  A"  and  B"  (in  the 
deformed  state)  is  therefore  given  by  the  vector  relation 


to  - tx  r 


(3.3b) 


s 


(3.4) 
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Fig.  3.3  Kinematics  of  surface  point  displacements  in 
conformal  contact. 
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Only  the  components  of  displacement  along  the  r axis  will 
be  considered  in  the  formulation  of  the  contact  criterion.  The 
elastic  displacements,  and  w2  form  angles  of  X and  tt-X  respectively 
with  the  r direction  where 


ir  - y - <f> 

Z 


(3.5) 


The  rigid  body  displacements  and  A2  form  angles  -a  and  ir-a 
respectively  with  the  r direction  while  f forms  an  angle  K with  the 

A 

r direction.  E,  is  determined  by  the  shape  of  the  contour  curve  at 

^ A 

points  A and  B.  The  projections  of  S in  the  r direction  may  be 
written  as 


Iflm.  f -6  JO04.<k  f-  ('•*/,  -tVVj  % 


(3.6) 


where  6 = -(A^  + A^) . 6 is  termed  the  approach  and  represents  the 

distance  that  points  on  one  body  move  along  the  z axis  towards 
points  on  the  other  body  due  to  rigid  body  movement. 

In  general  the  displacement  w-j  and  w2  may  be  written  in 
terms  of  the  interfacial  pressure  plx'.y'.z1)  and  an  influence 
function  G-|  (x,x' ,y,y ' ,z,z ' ) as 


J n. 


(3.7) 
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where  ft.  is  the  area  of  contact  of  body  i.  Physically  the  function 

.1 

G.j  ( x,x!y,y ' z,z' ) represents  the  elastic  displacement  in  the  -n  . 
direction  at  point  (x,y,z)  due  to  a unit  load  at  point  {x\y',z')  in 
the  -n^  direction.  In  the  most  general  form  equation  (3.6)  may  be 
written  as. 

St-  - If  I f " & 01  + 

+ XAi  A I Tf<«;  i)  Gf  (*,*',  iX  + 

4 /OH  aJ  n*w)  (3.8) 

In  order  to  solve  the  conformal  contact  problem  it  is 
necessary  to  find  the  interfacial  pressure  P(x',y',z'),  the  approach 
6,  and  the  final  contact  area  ft,  all  of  which  satisfy  equation  (3.8) 
and  the  following  boundary  conditions: 


S>  -o 

INSIDE  41 

(3.9a) 

TV*,?,?)  2 0 

INSIDE  SL 

(3.9b) 

S,  >0 

OUTSIDE  SL 

(3.9c) 

Pcx.i.i)  =0 

OUTSIDE  jCI 

(3.9d) 

Furthermore,  it  is  required  to  verify  the  accuracy  of  equation 
(3.2),  i.e.,  that  the  points  within  the  contact  area  satisfy  the 
relationships: 
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St  =o 

S«j  = o 


/first  p£  Jl 


/NS/Pe  SL 


(3.9e) 


(3.9f ) 


where  S.  and  S represent  the  separation  in  the  t and  w directions 

t 03 

respectively. 


criterion  for  conformal  contact.  They  are  analogous  to  equations 
(2.6)  and  (2.7a-d)  for  nonconformal  contact.  Equations  (3.8)  and 
(3.9a-d)  may  be  solved  in  a similar  manner  as  equations  (2.6)  and 
(2.7a-d)  using  the  "Simply-Discre  cized"  method  of  solution. 

However,  this  solution  only  guarantees  the  displacements  be  compat- 


ible in  the  r direction  since  only  Sr  was  involved  in  equation  (3.8). 
Therefore  it  is  necessary  to  examine  the  components  of  the  displace- 

A A 

ments  in  the  t and  w directions  to  insure  that  A and  B merge  as 
originally  assumed.  An  iterative  scheme,  termed  the  "point-mating 
procedure,"  is  outlined  which  shows  how  successive  "Simply-Discretized" 
solutions  may  be  utilized  to  converge  upon  a final  solution  in  which 
merging  points  on  each  body  have  been  identified. 

Consider  the  "point-mating  procedure"  on  the  first  attempt 
at  solution.  As  shown  in  figure  3.1 


where  k-j  = 1 . With  the  "Simply  Discretized"  method  of  solution, 
equation  (3.8)  is  written  for  N + 1 field  points,  thus  for  each 


Conditions  (3.8)  and(3.9a-f)  represent  the  contact 


? k,  x, 


(3.10) 
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field  point  equation  (3.10)  must  be  assumed.  This  can  be  restated  as 

* . 

4 1 - k,  a,‘  (i-i.Nti)  (3.ii) 

The  first  solution  of  equation  (3.8)  yields  a pressure  field, 
approach  and  a contact  region  ft.  The  elastic  displacements  w-j  and 
at  each  field  point  may  be  calculated  via  equation  (3.7).  Now 
it  remains  to  see  if  indeed  the  field  points  on  each  body  merge. 

In  order  to  check  the  final  separation  it  is  necessary  to  compute 

A /\ 

the  displacements  in  both  the  w and  t direction.  Denoting  the  dis- 
placements of  a field  point  on  body  i in  the  tangential  plane  of  body 
i by  u^  and  respectively,  they  may  be  determined  as  in  equations 
(3.12a ,b) , where  u^  lies  in  the  plane  of  the  contour  curve  and  v^  is 

A 

in  the  direction  of  w. 

Ui  ')  (*,*■',  'iA  ,1,1') dt'didl'  (3.12a) 

i i 1. iA W. if ' dl ' (3.12b) 

In  equations  (3.1 2a , b ) /-/■  and  I..  represent  the  influence  functions 
for  body  i which  relate  the  displacements  u . and  v.  respectively 
at  a point  (x,y,z)  to  the  normal  pressures  exerted  at  a point 
(x',y',z').  The  pressure  field  P(x',y',zl)  and  the  contact  areas 
ft.j  are  known  from  the  solution  of  equation  (3.8). 

Examining  one  set  of  points,  A on  body  1 and  B on  body  2, 
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which  are  assumed  to  merge  after  deformation,  the  total  separation 
may  be  written  as: 

S - St  f + St  * * Su  u = ) co*\  r + 

+(m,-ut)c<n\i  +(<f-if,)A>  *fxhfti -£rr  (3.13a) 

where 

f s fr  r * fu  to  (3.13b) 

and 

£ i = -(*,+£1)$  = iT  r +6ti  + O u>  (3.13c) 

Equation  (3.8),  which  was  originally  solved  by  the  "Simply-Discretized" 

A 

method  represents  the  r components  of  equation  (3.13). 

In  the  first  solution  of  equation  (3.8)  it  was  assumed  that 
u,  = u9,  v,  = v9,  f = 0.  However,  the  first  solution  has  now 

c.  \ C CO 

provided  values  for  u^ , u^,  v-j , and  6.  f is  a function  of  the 

geometry  of  the  undeformed  surfaces  and  6 may  be  decomposed  into  6r 
and  6^.  Therefore  all  the  quantities  in  equation  (3.13)  are  known 
to  a first  approximation  after  the  first  simply  discretized  solution. 
The  separation  after  the  first  solution  may  be  computed  and  in  general 
it  will  be  non-zero.  It  will  be  shown  how  the  separation  may  be 
utilized  to  better  approximate  equation  (3.11)  sc  that  the  separation 
of  a second  solution  to  the  same  problem  will  be  much  smaller. 
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In  order  to  simplify  the  illustration  of  the  point-mating 
procedure,  it  will  be  assumed  that  for  the  problem  at  hand  v-|  and 
are  zero  for  all  field  points.  This  is  the  case  for  axi symmetric 
bodies.  Furthermore,  consider  the  seat  to  be  fixed  at  some  point 
well  removed  from  the  contact  area  (where  the  elastic  deformation 
is  negligible),  thus  after  the  deformation  the  indentor  will  have 
displaced  the  entire  amount  6 due  to  rigid  body  displacement.  Typi- 
cal displacements  for  points  A and  B,  including  the  u-j  and  Ug 
displacements  found  via  the  first  simply  discretized  solution,  are 
illustrated  in  figure  3.4.  The  solution  to  equation  (3.8)  guaranteed 

A 

only  that  the  separation  in  the  r direction  would  be  zero.  As 
shown,  the  points  A and  B will  in  general  be  separated  by  a distance 
S^.  Equation  (3.11)  may  now  be  modified  such  that 

9 

% • 

- A}‘  A, * (3.14) 


(3.15) 


and  subscript  j refers  to  the  number  of  the  iteration. 

This  modification  compensates  for  the  error  in  the  original 
assumption  (3.11).  Using  relationship  (3.14),  the  calculation  may 
be  performed  a second  time.  The  value  of  S'  in  the  second  solution 
will  be  much  smaller  than  that  of  the  first,  however,  if  it  is 
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Fig.  3.4.  Displacements  of  typical  field  points  in 
conformal  contact. 
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still  too  large,  a new  value  of  k^,  may  be  computed  following 
equation  (3.15). 

Equations  (3.14)  and  (3.15)  may  be  generalized  as 


A 


where 


(3.16a) 


(3.16b) 


and  j denotes  values  associated  with  the  j th  "Simply  Discretized" 

solution  [k,  = 1].  The  value  of  k1.  in  equation  (3.16)  is  used 

J 

in  place  of  k1.  , for  the  j th  "Simply  Discretized"  solution.  The 
3 

iterative  scheme  can  be  repeated  until  the  desired  tolerances  on  Sj. 
are  met. 

A 

Also  of  concern  is  the  separation  in  the  u>  direction.  It 
was  assumed  that  points  along  contour  curves  merge,  however,  if  v-j 
and  ^2  are  n°t  equal,  points  A and  B may  not  lie  on  the  same  contour 
curve  after  deformation.  This  can  be  determined  by  examining  the 

separation.  S',  in  the  u>  direction.  It  must  be  remembered  that  the 

or 

solution  at  hand  is  an  iterative  one  and  therefore,  the  separation 
of  field  points  will  in  general  never  be  zero.  The  separation  can 
only  be  reduced  to  an  acceptable  amount. 

In  summary,  a mathematical  model  of  frictionless  conformal 
contact  theory  has  been  presented.  The  model  takes  the  form  of 


equations  (3.8)  and  (3.9a-d).  Since  there  is  no  knowledge  a priori 
of  which  points  on  each  surface  merge,  initially  points  located  by 
equation  (3.11)  are  assumed  to  merge.  The  contact  criterion  of 
equations  (3.8)  and  (3.9a-d)  only  insure  that  two  field  points, 
assumed  to  merge,  have  zero  separation  in  the  r direction.  In  order 
to  guarantee  the  absence  of  separation  between  two  field  points 
in  the  contact  region,  the  iterative  scheme  termed  the  "point-mating 
procedure"  must  be  applied. 


4.  GENERATION  OF  INFLUENCE  FUNCTIONS 


4.1  Introduction 

In  the  formulation  of  both  the  nonconformal  and  conformal 
contact  theories,  the  influence  function  plays  a crucial  role. 

Physically  the  influence  function  relates  the  elastic  dis- 
placement at  a given  point  to  the  applied  force  at  some  other  point. 
The  elastic  displacements  due  to  a given  pressure  distribution  can  be 
found  by  integrating  the  product  of  pressure  and  the  influence 
function  over  the  contact  area.  This  is  illustrated  by  equation  3.7. 

In  any  given  problem  it  is  necessary  to  know  the  influence 
functions  which  are  appropriate  for  the  given  surfaces  in  contact. 

In  nonconformal  theory  the  contact  area  is  approximated  by  a plane 
making  it  appropriate  to  use  the  Boussinesq  influence  function,  for  a 
point  load  on  a half  space,  as  the  influence  function  for  all  sur- 
face geometries.  However,  in  conformal  theory,  where  the  contact 
surface  can  not  be  approximated  by  a plane,  it  is  necessary  to  find 
the  influence  functions  explicitly  for  each  of  the  bodies  in  contact. 

For  some  problems  analytic  influence  funtions  may  be  found; 
however,  in  the  event  that  no  analytic  functions  are  available,  they 
may  be  generated  numerically.  The  following  two  sections,  4.2  and 
4.3,  present  the  classical  solutions  of  a half  space  loaded  under  a 
point  load  and  line  load  respecti vely . In  section  4.4  the  principles 


42 


43 


involved  in  the  numerical  generation  of  influence  functions  are 
developed.  The  remainder  of  chapter  4 contains  examples  of  influence 
functions  generated  with  finite  element  techniques.  Wherever 
feasible,  analytic  solutions  are  compared  to  the  numerical  influence 
functions. 


4.2  Influence  Function  for  a Point  Load  on  a Half  Space 

The  problem  of  a half  space  loaded  with  a normal  concentrated 
load  was  first  solved  by  J.  Boussinesq  [1885]  J Consider  the  half 
space  and  coordinate  system  illustrated  in  figure  4.1. 

A concentrated  load  F is  applied  at  point  0 and  point  A is 
located  on  the  surface  of  the  half  space  at  a distance  |r|  from  0. 
u represents  the  elastic  displacement  of  point  A in  the  direction  of 
r while  w represents  the  elastic  displacement  in  the  z direction. 
Boussinesq  found  u and  w to  be  given  by 


( l-2U)(  /■*■ u)  F 
Zrr  £ It  I 


(4.1) 


US 


F 

IT E ifl 


(4.2) 


where  E and  v are  Young's  modulus  and  Poisson's  ratio  respectively  for 
the  half  space. 

Equation  (4.2)  forms  the  basis  of  the  influence  function  used 
^See  Timoshenko  and  Goodier  [1970],  pp.  398-402. 
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Fig.  4.1.  Point  load  on  half  space 
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in  nonconformal  contact  theory.  We  wish  to  generate  a function 
G (x,  x',  y,  y')  which  relates  the  normal  displacement  at  (x,  y) 
to  a unit  load  at  (x',  y').  Noting  that  in  equation  (4.2)  w is  pro- 
portional to  l/|r|,  where  |rj  is  the  distance  between  points  (x,  y) 
and  (x',  y'),  the  influence  function  for  the  normal  surface  dis- 
placement on  a half  space  may  be  written  as 


(!-»')  i 

7T£ 


(4.3) 


Similarly,  for  the  displacements  tangent  to  the  surface,  a function 
H (x,  x',  y,  y')  may  be  defined  as 


~(l-2vXl+v)  / 

27 r£ 


(4.4) 


If  a pressure  field  p (x1,  y1)  were  considered  to  act  over 
the  surface  of  a half  space  within  a region  ft,  the  displacements  w 
and  u due  to  this  pressure  field  can  be  calculated  by  the  following 
equations : 

uj  --fTr*;y)6(xxir)4*'dt'  <4-5> 

%rt 

SI 


(4.6) 
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Both  G and  H play  an  important  role  in  the  generation  of  numerical 
influence  functions  which  will  be  explained  in  section  4.4. 


4.3  Analytic  Solution  for  a Line  Load  on  a Plane 

Consider  a line  load  acting  on  the  edge  of  a semi -infinitely 
plate.  It  is  desired  to  find  the  displacement  field  in  the  plate  due 
to  the  given  loading,  (see  fig.  4.2) 

As  posed  the  problem  is  one  of  plane  stress  and  was  originally 
solved  by  Flamant  [1892].  For  boundary  conditions,  it  is  assumed  that 
points  along  the  y axis  have  no  motion  in  the  x direction  while  a point 
A,  located  along  the  y axis  at  a distance  d from  the  surface,  is  fixed 
rigidly.  The  displacement  field  then  becomes^ 


(T-~ 


U s - 


2 P * (\-\>)F  _ . ^ 

— — — — — - Q AbnO  + 

7TB  V 7TB 

zvt 


(4.7) 
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(4.8) 


where  E and  v are  the  Young's  modulus  and  Poisson's  ratio  respectively 
of  the  plate.  The  solution  is  not  unique  in  that  it  is  dependent  on 
the  value  of  d chosen  in  the  boundary  condition. 
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See  Timoshenko  and  Goodier 


[1970],  p. 


103. 
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4 


Fig.  4.2.  Plate  loaded  under  line  load 
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Let  0 be  tt/2.  The  surface  displacements  for  the  plane  stress 
problem  become 


U s + 


ZF 

TT  £ 


O+v) 

TTf 


(4.9) 


l r - - 


a-v) 

ZE 


(4.10) 


For  the  case  of  plane  strain  equations  (4.9)  and  (4.10)  may  be  written 
substituting  v by  v and  E by  E where 


(4.11) 


The  proof  of  this  substitution  is  shown  in  appendix  1.  For  plane 
strain  the  surface  displacements  at  9 = tt/2  become 
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2F(l-»%)  . t 
7 TE  ^fd. 


(i+v) 


(4.12) 


_ (l-Z\t)(l+V)F 

ZE 


Equations  (4.12)  and  (4.13)  form  the  basis  of  an  influence 
function  for  a line  load  on  a plane.  Consider  a unit  line  load. 
Parallel  to  the  z axis  at  x',  then  the  displacements  along  a line  at 
coordinate  x becomes 


G(x.x') 


Z U-v*)  „ Ix-k'I 


U+v) 

ne 


(4.14) 


H (XX)  - - 


0-Z»)(\4U)  xX 
IE  lx-x'l 


(4.15) 


where  G (x,  x')  and  h(x,  x')  represent  the  influence  functions  for 
displacements  u and  v respectively  for  0 = tt/2.  Equations  (4.14)  and 
(4.15)  are  essential  to  the  numerical  generation  of  influence  functions 
when  line  loads  are  involved.  These  results  will  be  used  in  sections 
4.7  and  4.8  for  line  loads  on  bodies  with  cylindrical  surfaces. 
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4.4  Numerical  Generation  of  Influence  Functions 

Consider  a three  dimensional  body  whose  surface  is  defined 
by  z = <|>  (x,  y),  as  illustrated  in  figure  4.3.  For  a linear  isotropic 
material,  the  influence  function  for  the  surface  point  displacements 
will  be  of  the  form 


where  F is  a unit  load  acting  normal  to  the  surface  at  point 
(x‘,  y',  z1)  and  (x,  x',  y,  y1)  represents  the  displacement  in  a 
defined  direction  at  some  other  surface  point  (x,  y,  z).  E is  the 
Young's  modulus  of  the  material  and  v is  Poisson's  ratio.  The 
function  g(x,  x',  y,  y',  v)  depends  on  the  geometry  of  the  body  and 
has  dimensions  of  [1/L].  In  some  cases  g may  be  found  analytically 
however  for  more  complicated  geometries  the  task  may  be  impossible. 
When  g can  not  be  derived  by  analytic  means,  it  can  in  principle  be 
constructed  from  a set  of  finite  element  solutions.  A simple  example 
follows  which  demonstrates  this  numerical  procedure. 

Consider  the  region  (z  < 0)  defined  by  the  curvilinear 
coordinate  system  illustrated  in  figure  4.4. 

A unit  load  F is  applied  to  point  (x',  y')  normal  to  the 
surface  and  it  is  desired  to  find  the  displacement  Gina  specified 
direction  at  a point  (x,  y).  For  the  purpose  of  finding  finite 


(4.16) 


Q.C 


Fig.  4.3.  Three  dimensional  surface  4>  (x,  y)  under  point  load 


Fig.  4.4.  Curvilinear  coordinate  system  on  three  dimensional  surface 


52 


element  solutions,  a model  of  the  body  is  divided  into  "elements." 

The  surface  of  the  body  is  defined  by  the  top  sides  of  some  of  those 
elements.  The  intersection  of  more  than  two  elements  defines  a node. 
Figure  4.4  illustrates  a typical  discretization  of  the  surface  of 
some  typical  body  into  elements. 

In  order  to  find  G at  (x,  y)  due  to  F at  (x1 , y' ) it  will  be 
necessary  to  solve  a number  of  finite  element  problems  with  the  above 
model.  Consider  that  a set  of  solutions  is  known  via  finite  element 
techniques  within  the  region  of  interest  around  (x,  y)  and  (x',  y'). 
Each  solution  corresponds  to  a problem  where  the  point  load  is 
applied  to  a different  node.  Thus  from  these  solutions  the  displace- 
ment at  each  node  is  known  due  to  a point  load  at  any  of  the  other 
nodes. 

Now  consider  in  further  detail  elements  i and  j which 
contain  points  (x',  y1)  and  (x,  y)  respectively,  (see  fig.  4.4) 

Let  us  approximate  element  1,  2,  3,  4 as  a plane  facet,  then 
a set  of  four  forces  located  at  nodes  1,  2,  3 and  4 may  be  found 
which  is  equipollent  to  F at  (x‘,  y1)  i.e.,  if  it  is  required  that 
the  forces  at  nodes  1,  2,  3 and  4 sum  to  F and  that  their  moment  about 
(x1,  y')  is  zero  then  F may  be  replaced  by  forces  Fx , F2,  F3,  and  F4 
at  nodes  1,  2,  3,  4 respectively.  By  superposition  of  the  finite 
element  solutions,  the  displacement  at  node  5 may  be  found  due  to  the 
set  of  forces  Fx,  F2,  F3,  F4 . Similarly  the  displacement  at  nodes  6,  7, 
8 may  be  found.  One  final  interpolation  may  be  made  between  these  dis- 
placements at  nodes  5,  6,  7 and  8 to  obtain  the  displacement  G at  (x,  y) 
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due  to  F.  The  above  scheme,  involving  the  interpolation  between 
finite  element  solutions,  illustrates  how  an  influence  function  can 
be  generated  numerically  for  arbitrary  surfaces. 

It  must  be  noted  that  finite  element  solutions  for  a point 
load  applied  to  a node  will  yield  a finite  displacement  directly 
under  the  load.  This  is  inconsistent  with,  at  least,  equations  (4.1) 
and  (4.2)  which  predict  infinite  displacements  under  a point  load  on 
a plane.  The  finite  element  displacement  function  is  only  valid  away 
from  the  point  load.  In  the  neighborhood  of  the  load,  the  appropriate 
singularities  must  be  identified  as  described  in  sections  4.5--4.8. 

The  following  sections  contain  examples  of  the  generation  of 
influence  functions  via  the  method  described  above.  They  deal  with 
both  point  and  line  loading  on  spherical  and  cylindrical  surfaces 
respectively.  In  all  but  one  of  the  examples  analytic  solutions  are 
compared  to  the  numerically  generated  influence  functions. 

4.5  Influence  Function  for  a Point  Load  on  a Sphere 

Sternberg  and  Rosenthal  [19521  have  found  the  solution  for 
the  stress  distribution  in  an  elastic  sphere  under  two  equal  and 
diametrically  opposed  point  loads.  Guerrero  and  Turteltaub  [1972] 
have  analyzed  a similar  problem  of  a sphere  under  a finite  number  of 
concentrated  surface  loads  of  arbitrary  orientation.  Both  of  these 
solutions  are  useful  in  providing  an  analytic  influence  function  for 
a point  load  on  a sphere.  The  results  that  follow  are  from  the 
analysis  of  Sternberg  and  Rosenthal. 

Consider  a sphere  compressed  by  two  concentrated  forces  F as 
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shown  in  figure  4.5.  From  dimensional  considerations,  symmetry 
considerations  and  the  fact  that  the  displacements  must  be  proportion- 
al to  loads  in  classical  problems  of  elasticity,  it  follows  that 
displacements  on  the  surface  of  a sphere  must  be  of  the  form 

T v> 


It  is  shown  in  Lure'  [1964]  that  the  displacements  u (radial)  and  uQ 

r u 

(meridional)  on  the  surface  of  the  sphere  are  as  follows: 


(m~%)  F F f iw-i  . ! @ 
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Fig.  4.5.  Sphere  under  diametrically  opposed  point  loads 


56 


uu 


f HI 


■ 


@ . B 

t f 2 


«*44m  © 


» m2+8in-3 

* — — — K 


m 


-/«<*if  )mt9+  ^ M*%  dJUtij  *&(  j —— — i— j 


, £ o dPm(M@) 

+ <£  

*?/ 


(4.19) 


where  m is  the  reciprocal  of  Poisson's  ratio,  G is  the  modulus  of 
rigidity  and  (cos  6)  are  the  Legendre  polynomials  in  cos  0.  The 
coefficients  are  given  as 
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m 3'(n-i)(n  +0(n+ 2 j A' 


[(7ff\3  + ll  m2  -28m  ■tii>)(h+2)2’ 

+ (-I8m3- 46/n*  + l3Zm -72Xn*2) 

+ y(m-i)(mz  +8m -8 ) ] 


where 


(4.21) 


nfn-/)  +(in  + i) 


m+i 

m 


(4.22) 


The  influence  functions  for  a point  load  on  a sphere  may  be 
obtained  by  considering  a unit  load  F in  equations  (4.18,  4.19)  and 
letting  9 represent  the  angle  between  the  vectors  describing  the 
position  of  the  load  and  the  point  A where  the  displacements  are 
desired.  This  is  illustrated  in  figure  4.6. 

Consider  points  A and  B on  the  surface  of  a sphere.  A 

AAA 

cartesian  coordinate  system  (n,  y>)  is  constructed  with  the  origin, 

0,  at  the  center  of  the  sphere.  The  vectors  0A  and  OB  form  an  angle  9 
between  them. 

If  loads  F are  applied  at  point  B and  at  a point  diametrical- 
ly opposite  D then  the  displacements  u and  u.  at  A are  defined  by 

r 0 

equations  (4.18)  and  (4.19)  respectively.  It  is  important  to  note 
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Fig.  4.6.  Coordinate  system  for  influence  functions  for  a 
point  load  on  a sphere. 
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that  the  direction  of  uQ  is  not  only  a function  of  6 but  also  of  the 
spherical  coordinates  (iK  6)  of  D and  those  (<K , S')  of  A.  uQ  is 
measured  along  the  tangent  to  the  great  circle  at  D which  passes 
through  point  A. 

The  influence  functions  for  a point  load  on  a sphere  may  be 
written  as 


<sc% fi, y'J'a £,*)--  (*>») 


(4.23) 


and 


(4.24) 


where 


&(V, £ y'/s')  =c^e/xn. 


(4.25) 
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| OD 


From  spherical  law  of  cosines  or  from  cos  0 = OA  • 0D/|0A|* 
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G and  H represent  the  influence  functions  for  u^  and  u0  respectively 
which  are  given  by  equations  (4.18)  and  (4.19)  respectively.  The 
displacements  u^  (ip,  6)  and  uQ  (iJj,  6)  at  point  (^,  6)  due  to  a dis- 
tributed load  p (ip ' , 8')  over  region  ft  can  be  calculated  respectively 

4*ta/0  -■/  *W)  octjM »,*,* 

A 


and 


ctgfyii):]*  t(ry)H(rpyjy,eJii)R,Ainr'd<f'dr'  <4.27) 


Now  consider  the  numerical  generation  of  these  influence 
functions.  It  is  known  from  symmetry  that  the  magnitudes  of  G and  H 
are  related  to  R,  0,  v and  E.  Furthermore  the  orientation  of  positive 
G is  always  radially  inward  whereas  that  of  H is  tangent  to  the  sur- 
face in  the  direction  away  from  the  unit  load.  Therefore  it  is 
necessary  to  find  functions  gr  (e,  v)  and  g0  (9,  v)  in  equations 
(4.23)  and  (4.24).  These  functions  can  be  constructed  if  desired^ 
from  one  finite  element  solution  of  a point  load  on  a sphere. 


The  analytic  solution  of  Sternberg  and  Rosenthal  [1952]  will 
serve  as  a check  on  the  finite  element  development  which  will  then  be 
applied  to  the  spherical  seat  and  other  problems  without  analytic 
solutions. 
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Consider  the  sphere  under  the  loading  in  figure  4.5.  This 
problem  can  be  modeled  for  finite  element  analysis  with  the  grid 
shown  in  figure  4.7.  The  present  analysis  was  made  with  the  computer 
program  for  axisymmetric  problems,  described  in  Wilson  [19651.  The 
grid  in  figure  4.7  represents  a cross  section  of  ring  elements  which 
are  axisymmetric  about  the  z axis.  Because  of  the  symmetry  of  the 
loading,  points  located  on  the  x axis  were  restricted  to  move  only 
in  the  x direction  while  points  on  the  z axis  were  fixed  from  moving 
in  the  x direction.  The  output  data  of  interest  are  the  displacements 
of  the  surface  points.  Consider  the  finite  element  analysis  of  a 

C 

sphere  where  E = 30  x 10  psi,  v = 0.3,  R = 1 in.  which  is  compressed 
between  two  forces  F,  where  F = 30  x 10^  lb.  From  equation  (4.17)  it 
follows  that 


and  the  function  h and  hQ  are  exactly  those,  g and  gQ,  in  equations 

r u r u 

(4.23)  and  (4.24). 

The  functions  h (0,  v)  and  hQ  (0,  v)  are  non-dimensional 
r u 

functions  of  displacement  which  are  known  at  the  nodes  in  figure  4.7 
from  the  finite  element  analysis.  Knowing  the  displacements  at  the 


(4.28) 


(4.29) 
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A 


Fig.  4.7.  Axi symmetric  finite  element  model  of  a sphere  under 
two  diametrically  opposed  point  loads. 


Fig.  4.8.  Matching  of  analytic  singularity  to  numerical 
(finite  element)  influence  function. 


COMPARISON  OF  INFLUENCE  FUNCTIONS  GENERATED  USING  FINITE  ELEMENTS 
TO  THE  ANALYTIC  SOLUTIONS  FOR  POINT  LOADS  ON  A SPHERE 
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surface  nodes,  the  displacement  at  any  surface  point  may  be  obtainea 
through  interpolation. 

The  computed  displacements  under  the  load  are  finite  and  the 
numerical  results  must  be  matched  at  some  point  near  the  load  with  the 
singularity  of  the  function.  It  is  shown  in  appendix  A that  the 
dominant  singularities  in  the  displacements  under  a point  load  on  a 
sphere  are  indeed  those  of  equations  (4.1)  and  (4.2). 

It  was  found  that  with  an  appropriately  dense  mesh,  such  as 
that  shown  in  figure  4.7,  the  numerical  influence  function  hr  (0,  v) 
would  merge  with  its  analytic  singularity  over  a region  near  the  load. 
This  is  illustrated  in  figure  4.8  where  the  numerical  influence  func- 
tion hr  merges  with  the  singular  function  in  the  region  between  A and 
B.  The  singularity  was  matched  to  the  numerical  function  at  A,  thus 
between  0 and  0^  the  behavior  of  the  generated  influence  function  was 
taken  to  be  that  of  the  singularity  while  for  0 greater  than  0^  the 
numerical  values  of  hr  were  utilized  to  describe  the  influence 
function  gr  (0,  v).  Similar  treatment  was  used  to  generate 
g0  (e,  v). 

A comparison  of  the  displacement  functions  given  by  "analytic" 
equations  (4.18)  and  (4.19)  with  the  displacement  function  hp  (0,  v) 
at  the  node  points  of  finite  element  analysis  are  given  in  table  4.1. 

It  can  be  concluded  that  the  numerical  influence  functions  are  accurate 
representations  of  the  analytic  ones. 

4. 6 Influence  Function  for  a Point  Load  on  a Spherical  Cavity 

The  displacements  due  to  a point  load  on  a spherical  cavity 
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have  not  been  published  to  date.  An  investigation  by  Sternberg  et  al 
[1951],  which  deals  with  the  solution  to  the  axisymmetric  problem  of 
a region  bounded  by  two  concentric  spheres,  could  be  used  under 
certain  limiting  conditions  to  produce  the  influence  function;  how- 
ever, the  limit  process  is  very  involved  and  has  not  been  performed. 
For  the  purposes  of  the  present  research,  the  influence  function  for 
a point  load  on  a spherical  cavity  is  derived  numerically  as  outlined 
in  this  section. 

Consider  a spherical  cavity  of  Radius  R under  two  diametri- 
cally opposed  point  loads,  F.  The  material  has  elastic  moduli  E and 
v.  (see  fig.  4.9) 

The  displacements  ur  and  uQ  are  desired  as  a function  of  R, 

0,  E,  v,  and  F.  From  equation  (4.17)  which  is  also  valid  for  a 
spherical  seat,  we  note  that  the  displacements  must  be  at  the  form 


(4.30a) 


and 


(4.30b) 
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Fig.  4.9.  Spherical  cavity  under  diametrically  opposed  point  loads 
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Now  consider  a (r^,  T2)  coordinate  system  as  shown  in  figure  4.11. 

The  influence  function  for  radial  and  tangential  displace- 
ments respectively  due  to  loads,  diametrically  opposed  on  a spherical 
cavity,  may  be  expressed  in  terms  of  a ^ » T2)  sPherica^ 
coordinate  system  as 


We  seek  to  find  the  functions  gr  and  gQ  which  are  the 
displacements  on  the  internal  spherical  boundary  surface  of  an  infinite 
region  under  a loading  condition  such  that  F/ER  = 1.  These  displace- 
ment functions  can  be  easily  found  using  a single  finite  element 
analysis. 

Consider  the  discretized  model  of  a spherical  cavity  in  an 
infinite  medium  as  shown  in  figure  4.11.  The  model  represents  the 
loading  illustrated  in  figure  4.9.  Each  element  is  a ring, 
axisymmetric  about  the  z axis.  The  boundary  conditions  for  the  model 
restrict  the  nodes  on  the  z axis  to  move  only  in  the  z direction  while 
those  on  the  x axis  are  allowed  to  move  only  in  the  x direction.  The 


(4.31) 


and 


% (6,  v) 


(4.32) 
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Fig.  4.10.  Coordinate  system  for  influence  functions  for  a 
point  load  on  a spherical  seat. 
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nodes  on  the  outer  boundary  were  free.  The  values  of  the  Force  F, 
radii,  Poisson's  Ratio  and  Young's  modulus  used  in  the  finite  element 
analysis  were  as  follows 

F = 30  x io7  lb. 

R = 10  in. 

R = 200  in. 
o 

v = 0.3  and  0.25 

E = 30  x io6  psi  (4.33) 

The  solution  of  the  finite  element  analysis  gives  the  values  of  the 

displacements  on  the  surface  of  the  cavity  at  the  nodal  points  in 

figure  4.11.  These  displacements,  u^  and  uQ,  under  conditions  where 

F/ER  = 1,  actually  represent  the  values  of  the  functions  gr  (9,  v)  and 

g0  (0,  v).  Therefore  by  interpolating  between  these  values  an 

approximation  of  g^  (6,  v)  and  g_  (9,  v)  is  known  for  all  0.  As  with 

r w 

the  influence  functions  for  the  sphere,  the  singularities  of  equations 
(4.31)  and  (4.32)  near  the  point  load  were  represented  by  equations 
(4.1)  and  (4.2). 

4.7  Influence  Functions  for  a Cylinder  Under  Concentrated  Line  Loads 
Consider  a long  cylinder  under  two  concentrated  line  loads  as 
shown  in  figure  4.12. 

The  problem  is  one  of  plane  strain  and  has  been  solved  in 
Muskhel ishvil i [1963].  The  displacements  at  the  point  Q,  in  terms  of 
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quantities  labeled  in  figure  4.12,  are 


and 


(4.34) 


- (jUn  2c(,  + Ao» i 


(Mm)R  3 . 


(4.35) 


where  F is  a line  load  expressed  in  units  of  force/length,  ux  and  u^ 


are  the  displacements  in  the  x and  y directions  respectively,  andAand 
y are  the  Lame'  constants  of  the  cylinder.  Equations  (4.34)  and 
(4.35)  may  be  combined  to  find  the  displacement  field  in  polar 
cordi nates  of  a cylinder  under  two  diametrically  opposed  line  loads, 
as  illustrated  in  figure  4.13. 

The  derivation  of  ur  and  uQ  is  performed  in  appendix  B. 

The  results  are  as  follows: 


Uy  - - f ^ j^/Co4  OJIm.  /ah  * / j ir  F 


JjT  Mm6  (4.36) 


Ak 
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A 


Fig.  4.13.  Polar  coordinate  system  for  a cylinder  under  two 
diametrically  opposed  line  loads. 


and 
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Ue  - - F X,  M*<-8  Un  | - P ~f  (4.37) 


where 


2(1- V1) 

ir£ 


(4.38) 


and 


(i+v)(i-zv) 

2£ 


(4.39) 


The  influence  function  for  a line  load  on  a cylinder  may  be  formed 
from  equations  (4.36)  and  (4.37)  by  considering  0 to  be  the  angle 
between  the  vector  describing  the  position  of  the  unit  load  F and  the 
point  Q.  The  influence  functions  in  the  radial  and  tangential 
directions  respectively  are 


G ( Gj  Vj  £ ) r “ Kf  (B  -B  )£vl  fan  f / J 

Kz 


f ~jr  (0- 0 ) 


(4.40) 
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and 


H( 6)6>  Vj  E)  <4**  (&-6')jUt  fa*.  [-g^] 

+*  if  Kl  /C*l(&-Q')  (4.41) 


where  6 and  0'  are  illustrated  in  figure  4.14. 

It  is  important  to  note  that  equations  (4.40)  and  (4.44) 
represent  the  influence  function  for  two  diametrically  opposed  line 
loads.  They  may  be  integrated  as  shown  below  to  find  the  displace- 
ments ur  (0,  R)  and  uQ  (0,  R)  due  to  two  symmetric  loadings  p (0), 
where  p (0)  = p (0  + tt).  The  displacements  may  be  calculated  by 

UrOj  £)  p(e  ')Rde'  (4.4 

and 


uecs>/?)*f  ej  v,  e)  <f(e ')  Rde ' 

A 


(4.43) 


Now  consider  the  numerical  generation  of  equations  (4.40)  and 
(4.41).  It  is  known  that  the  displacements  up  and  u0  due  to  diametri- 
cally opposed  point  loads,  F are  functions  of  0,  F,  E,  v and  R,  where 
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Fig.  4.14.  Coordinate  system  for  the  influence  functions  of 
a cylinder  under  two  diametrically  opposed  line  loads. 


I 
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E and  v are  the  elastic  constants  of  the  cylinder  of  radius  R.  From 

dimensional  analysis  the  displacements  u and  u„  may  be  expressed  in 

r y 

the  forms 


(4.44) 


T * *'■£'**’> 


Furthermore,  because  of  the  linearity  of  the  problem,  the 
ratio  er  must  appear  linearly  in  the  above  equations;  hence  ur  and 
uQ  may  be  expressed  in  the  form: 


(4.46) 


and 


(4.47) 


We  therefore  seek  to  find  the  dimensionless  functions  g (0,  v)  and 
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g0  (0,  v)  which  represent  the  displacements  u^  and  uQ  under  the 
conditions  that  F/E  = 1.  These  functions  may  be  found  with  one  finite 
element  analysis  of  a cylinder  under  two  diametrically  opposed  line 
loads. 

Consider  the  model  of  a cylinder  given  by  the  finite  element 
representation  in  figure  4.15.  The  elements  are  plane  stress  elements. 
The  line  load  per  unit  thickness,  F,  is  applied  and  due  to  the 
symmetry  of  loading  only  the  upper-right  quarter  of  the  cylinder  is 
considered.  Nodes  on  the  x axis  are  restricted  to  move  only  along 
the  x axis  and  nodes  on  the  z axis  are  restricted  to  move  only  along 
the  z axis.  For  the  solution  at  hand  the  values  of  the  parameters 
P,  E,  v and  R were  taken  as 

P = 3 x 107  lb/in. 

E = 32.967  x IQ6  psi 

(4.48) 

v = 0.42857 

R = 1 in. 


The  values  of  E and  v used  above  represent  the  moduli  in  plane  stress 
which  are  equivalent  to  the  values  of  E = 30  * 106  psi  and  v = 0.3  in 
plane  strain.  These  values  were  determined  via  equation  (4.11)  and 
the  results  are  appropriate  for  the  plane  strain  model  with  the  noted 
values  of  E and  v.  The  solution  yields  the  values  of  the  displace- 
ments, ur  and  u0  at  the  surface  nodes  in  figure  4.15.  These 
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A 


Fig.  4.15.  Plane  stress  finite  element  model  of  a cylinder 
under  line  loads. 


TABLE  4.2 

COMPARISON  OF  INFLUENCE  FUNCTIONS  GENERATED  USING  FINITE  ELEMENTS 
TO  THE  ANALYTIC  SOLUTIONS  FOR  LINE  LOADS  ON  A CYLINDER 
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displacements  represent  the  values  of  gr  (0,  v)  and  gQ  (0,  v)  at  the 
nodal  points  in  the  finite  element  model.  The  functions  gr  (0,  v) 

and  g (0,  v)  can  therefore  be  approximated  by  interpolating  between 

0 

these  values. 

The  singularities  for  the  functions  gr  and  gQ  are  given  by 
equations  (4.14)  and  (4.15),  that  is,  those  singularities  appropriate 
to  the  line  load  on  a plane.  It  is  shown  in  Appendix  C that  these 
are  indeed  the  correct  singularities  since  the  functions  in  equations 
(4.12)  and  (4.13)  approach  infinity  in  exactly  the  same  manner  as 
equations  (4.36)  and  (4.37),  respectively. 

It  may  be  concluded  that  by  knowing  values  of  gr  (0,  v)  and 
g0  (0,  v)  at  the  nodal  points  in  figure  4.15  and  by  having  correctly 
identified  the  singularities  of  those  functions,  the  influence 
functions  for  line  loads  on  a cylinder  have  been  determined.  As  a 
final  check  on  the  accuracy  of  such  a function  the  values  of  the 
displacements  g^,  and  g0  from  the  finite  element  analysis  have  been 
compared  to  those  of  equations  (4.12)  and  (4.13)  in  table  4.2.  It  is 
clearly  indicated  that  the  finite  element  solution  gives  a very 
accurate  representation  of  the  displacements  up  to  within  a half 
degree  from  the  applied  force. 

4.8  Influence  Function  for  a Line  Load  on  a Cylindrical  Cavity 

Consider  a cylindrical  cavity  within  an  infinite  solid  body. 
Given  that  two  line  loads,  F,  diametrically  opposed,  are  applied  to 
the  cavity  surface,  the  displacements  ur  and  u0  are  derived  from  known 
elasticity  solutions  in  Appendix  D.  (see  figure  4.16) 


Fig.  4.16.  Cylindrical  cavity  under  two  diametrical ly 
opposed  line  loads. 
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They  are  as  follows: 


Ur  r ft K,  <440  jbtfan 


(4.49) 


and 


Lj#  r F £ AJ  &m,  2.  ~ <04*  oj  (4.50) 


where 


ad-»z) 

7 T £ 


(4.51) 


and 


(l-  2.v)(l+  u) 

17 


(4.52) 


The  influence  functions  for  ur  and  uQ  can  be  derived  simply  from 
equations  (4.49)  and  (4.50).  By  considering  F to  be  a unit  load  per 
unit  width  and  by  replacing  0 by  9-9'  as  illustrated  in  figure  4.17. 


Fig.  4.17.  Coordinate  system  for  the  influence  functions  of 
a cylindrical  cavity  under  two  diametrically  opposed  line  loads. 
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the  influence  function  for  two  line  loads  on  a cylinder  become 


€(&,&', VtE)  --  -K,  (4.53) 

- Kt  M*.l(e-&)l 


It  is  interesting  to  note  that  the  influence  functions  for 
both  the  cylinder  and  cylindrical  seat  are  independent  of  the  radius 
R,  since  the  dimensional  analysis  leading  to  equations  (4.44)  and 
(4.45)  are  valid  for  the  cylindrical  cavity  as  well  as  for  the  solid 
cylinder.  Accordingly, the  displacement  functions  for  the  cylindrical 
seat  are: 


and 


H(0, B',U,E)  - I 


(4.54) 


(4.55) 


and 


(4.56) 
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where  $ represents  the  difference  between  the  angular  position  of  the 
force  and  the  displacement,  i.e.,  9-6'  in  figure  4.17.  The  functions 
9p  (<J>»  v)  and  gQ  (4>,  v)  represent  the  displacement  on  a cylinder  where 
the  ratio  F/e  =1.  An  approximation  to  these  functions  may  be  obtained 
through  finite  element  analysis.  From  equations  (4.55)  and  (4.56)  we 
see  that  the  required  influence  functions  are  found  by  dividing  g^  and 
9*  bV  E- 

Consider  the  model  of  a cylindrical  cavity  as  illustrated  in 
figure  4.18.  The  model  represents  the  loadings  illustrated  in  figure 
4.17.  The  points  along  the  z axis  are  restricted  to  move  only  along 
the  z axis  and  1 i kewise  the  points  on  the  x axis  are  restricted  to  move 
on  the  x axis.  The  outside  radius  is  thought  to  be  sufficiently 
removed  from  the  cavity  to  consider  it  at  infinity.  In  the  mathemati- 
cal problem  both  the  stresses  and  displacements  vanish  at  infinity. 
Therefore,  two  cases  of  boundary  conditions  will  be  considered  at  the 
outside  radius  in  the  discretized  model.  In  the  first  case  those 
points  are  free  while  in  the  second  they  will  be  fixed  rigidly. 

For  the  example  at  hand  the  values  of  E,  v,  R,  Rq  and  F were 
as  follows: 

F = 3 x io7  lb 
E = 32.967  x 106  psi 
v = 0.42857 
R = 1 in 
Ro  = 20  in 


(4.57) 


PLANE  OF  SYMMETRY 
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The  results  of  the  surface  nodal  point  displacement  for  the 
free  and  fixed  boundary  conditions  are  compared  to  the  analytic 
solutions  in  tables  4.3  and  4.4  respectively.  The  results  indicate 
that  the  boundary  condition  at  the  outside  radius  which  most  closely 
models  the  true  condition  is  when  those  points  are  considered  to  be 
free.  The  finite  element  data  in  table  4.3  shows  good  agreement 
between  the  numerical  and  analytic  solutions.  While  the  data  in 
table  4.4  also  indicates  a close  correspondence  it  does  not  agree  as 
well  as  the  results  in  table  4.3. 

Again  there  is  the  problem  of  finding  the  appropriate 
singularity.  As  might  be  expected,  the  Flamant  solution  for  a line 
load  on  a plane,  represents  the  singularity  for  the  case  of  a line 
load  on  a cylindrical  cavity.  This  is  proven  in  appendix  E. 

From  the  above  finite  element  analysis,  the  functions 
gr  (<J>,  v)  and  g^  (<t>,  v)  may  be  evaluated  away  from  the  applied  loads 
by  interpolating  between  the  nodal  displacements.  The  Flamant 
singularity,  given  by  equation  (4.14),  may  be  used  as  the  singularity 
of  the  function  (4.55),  while  the  constant  in  equation  (4.15)  may  be 
substituted  for  the  limit  of  equation  (4.56)  near  the  applied  loads. 
Having  constructed  equations  (4.55)  and  (4.56),  the  displacements  ur 
and  Uq  due  to  distributed  line  loads  p (0)  and  p (0  + it),  where 
p (0)  = p (6  + u),  may  be  calculated  from  equations  (4.42)  and  (4.43) 
respectively. 


TABLE  4.3 

COMPARISON  OF  INFLUENCE  FUNCTIONS  GENERATED  USING  FINITE  ELEMENTS  TO  THE 
ANALYTIC  SOLUTIONS  FOR  LINE  LOADS  ON  A CYLINDRICAL  CAVITY 
(OUTER  BOUNDARY  OF  FINITE  ELEMENT  MODEL  FREE) 
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COMPARISON  OF  INFLUENCE  FUNCTIONS  GENERATED  USING  FINITE  ELEMENTS  TO  THE 
ANALYTIC  SOLUTIONS  FOR  LINE  LOADS  ON  A CYLINDRICAL  CAVITY 
(OUTER  BOUNDARY  OF  FINITE  ELEMENT  MODEL  FIXED) 
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5.  CONFORMAL  ELASTIC  CONTACT  OF  A SPHERE 
INDENTING  A SPHERICAL  CAVITY 


5.1  Introduction 

Hertz  [1881]  provides  an  analytic  solution  to  the  problem  of 
two  contacting  bodies  with  quadratic  surfaces.  In  his  theory.  Hertz 
assumed  non-conformal  contact,  i.e.,  that  all  of  the  dimensions  of 
the  contact  patch  were  small  compared  to  the  radii  of  curvature  of  the 
bodies,  hence  justifying  the  approximation  of  the  surfaces  in  the 
contact  region,  by  elastic  half  spaces.  Hertzian  theory  may  be 
applied  to  the  problem  of  a sphere  indenting  a spherical  cavity,  how- 
ever, it  is  restricted  in  its  application  to  the  analysis  of  cases 
where  the  contact  patch  remains  small.  If  the  sphere  and  seat  are 
closely  conforming,  i.e.,  their  radii  are  nearly  equal,  then  Hertzian 
analysis  can  only  be  applied  for  very  small  loads  for  which  the 
contact  patch  remains  small. 

Goodman  and  Keer  [1965]  have  analyzed  the  conformal  problem 
of  a sphere  and  seat  with  nearly  equal  radii.  They  present  results 
for  areas  of  contact  larger  than  those  that  could  be  analyzed  by 
Hertzian  theory.  As  in  the  theory  of  Hertz,  Goodman  and  Keer 
assumed  that  points  on  surfaces  of  the  sphere  and  seat  respectively 
which  lie  along  a line  parallel  to  the  line  of  applied  load  will 
merge  after  deformation.  Furthermore,  the  basic  equation  in  their 
formulation  enforces  this  constraint.  It  will  be  shown  that  this 
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assumption  in  conformal  theory  can  lead  to  erroneous  strains  and  dis- 
placements when  large  areas  of  contact  are  analyzed. 

In  section  5.2  the  problem  of  an  elastic  sphere  indenting  an 
elastic  spherical  seat  will  be  formulated,  using  the  conformal  contact 
model  discussed  in  chapter  3.  The  numerical  procedures  used  in  the 
solution  will  be  reviewed  in  section  5.3  and  the  numerical  results 
will  be  presented  in  section  5.4. 

5.2  Formulation 

Consider  the  conformal  contact  problem  of  an  elastic  sphere 
of  radius  R-j  indenting  an  elastic  spherical  seat  of  radius 
(R2  > R-|).  It  is  assumed  that  the  sphere  and  seat  are  equilibrated  by 
pressure  distributions  equal  to  the  interfacial  contact  pressure  and 
diametrically  opposite  the  contact  region.  Contact  regions  will  be 
therefore  limited  to  hemispherical  contact.  A cross  section  of  a 
sphere,  body  1,  and  spherical  seat,  body  2,  in  point  contact  at  0 is 
illustrated  in  figure  5.1. 

AAA 

Let  us  establish  a global  coordinate  system  (x,  y,  z)  whose 
origin  will  be  fixed  at  the  initial  point  of  contact,  0,  such  that  the 
x - y plane  is  tangent  to  the  sphere  at  0 and  z is  directed  into  the 
sphere.  The  cross  section  of  the  sphere  and  seat  in  the  x - z plane 
represent  the  "contour  curves"  of  the  sphere  and  seat  and  each  curve 
is  a circle. 

It  is  initially  assumed  that  two  points  which  will  merge 
after  deformation,  A on  body  1 and  B on  body  2,  are  defined  such  that 
their  distances  from  the  origin  along  their  respective  contour  curves 
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Fig.  5.1.  Conformal  contact  between  a sphere  and  spherical 
seat  or  a cylinder  and  a cylindrical  seat. 
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are  equal,  i.e.,  = S2  (Fig.  5.1).  Points  A and  B are  defined  by 

vectors  r^  and  and  the  vector  difference  between  them  is  defined 

by  f,  the  vector  valued  profile  function.  Alternatively,  points  A 
and  B may  be  located  on  their  respective  contour  curves  by  angles 
measured  at  the  center  of  the  contour  curves  from  the  z axis  as 
shown  in  figure  5.1.  Point  A is  defined  by  angle  ip  and  point  B is 
defined  by  angle  4>.  In  the  notation  to  follow  in  this  chapter  all 
angles  \p  will  refer  to  locations  of  points  on  the  sphere  while  all 
angles  will  denote  points  on  the  seat. 

/\  A A 

A local  coordinate  system  (r,  t,  u>)  will  be  constructed  at 
point  A such  that  unit  vector  r,  which  represents  the  "mean  normal" 
forms  an  acute  angle  a with  the  z axis,  where 

e*  » -X - (5.i) 

£ 


unit  vector  t is  defined  to  be  tt/2  radians  clockwise  of  r and  w is 
such  that 


CO 


+ A 

t x r 


(5.2) 


The  contact  criterion  can  then  be  formulated  by  examining  the 
displacements  of  points  A and  B.  Consider  the  following  displacements 
as  illustrated  in  figure  5.2.  Point  A moves  through  a rigid  body 
translation  A^,  parallel  to  the  z axis, to  A1,  and  through  an  elastic 
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A* 


Fig.  5.2.  Displacements  on  sphere  and  spherical  seat 
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displacement  along  the  inward  normal  of  the  sphere  to  A". 

Similarly  point  B undergoes  a rigid  body  translation  parallel  to 
the  z axis  to  B1  and  an  elastic  displacement  w^  directed  along  the 
inward  normal  of  the  seat  to  B".  Because  of  the  symmetry  of  the 
contacting  surfaces  about  the  z axis,  there  will  be  no  displacements 
on  either  body  in  the  to  direction. 

We  now  impose  the  contact  criterion,  i.e.,  that  the 
projection  of  the  separation  (S  = A"  B")  in  the  direction  of  the  mean 

A 

normal  r must  vanish.  Following  equation  (3.6)  this  criterion  may  be 
expressed  in  the  form: 

Sf  - f 404  ? - £ 404  W*  )<C04.  \ (5-3) 


where 


and 


z 


S = - (a,  +a2 ) 


OAct&sn. 


/ fix  ~ \ 

' r,i  - rit  I 


(5.4) 

(5.5) 


(5.6) 


£ is  derived  in  appendix  F and  r 
components  of  vectors  r-j  and  r^ 


lx’  r2x’  rlz’ 
as  shown  in  fi 


r2z  are  x an<^  1 
gure  5.2.  f and  6 are 


the  profile  function  and  approach  respectively. 
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It  should  be  noted  that  the  peak  strains  will  remain  small 
(as  required  by  the  assumption  of  linear  elasticity  theory)  only  if 
the  initial  separation  f (see  fig.  5.1)  is  small  compared  to  the  local 
radii  of  curvature.  Therefore  some  additional  but  consistent  approx- 
imations are  possible  in  the  derivation  of  equation  (5.3),  (see 
appendix  H).  However,  these  approximations  result  in  no  significant 
reduction  of  computational  effort. 

Consider  now  the  (n-j » C-j  > y-j ) coordinate  fixed  to  the  sphere 
as  shown  in  figure  4.7.  Point  A can  be  located  in  spherical  coordi- 
nates by  (ty,  8,  R-j).  The  elastic  displacement  w-j  at  A,  shown  as  ur 
in  figure  4.7,  can  be  expressed  in  terms  of  the  pressure  distribution 
p (ip 1 , 8‘ ) in  the  form 

EJAnrdfi'dr'  (5-7) 
A 

where  G (ijj,  6,  xp ' , 8',  v-|,  E-j)  is  given  by  G in  equation  (4.22)  and 
n,  is  the  contact  surface  defined  on  the  sphere. 

Similarly  in  terms  of  the  (ri2>  coorcllnate  system  in 

figure  4.11,  W2  can  be  expressed  in  terms  of  p (<f>,  8)  by 

a 


(5.8) 
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where  (<j>,  3,  <J>',  3',  E^)  is  defined  by  equation  (4.35)  and  fig 

is  the  contact  are  defined  on  the  spherical  seat.  The  integrands 
in  equations  (5.7)  and  (5.8)  may  be  simplified  by  utilizing  the 
symmetry  of  the  pressure  field.  Because  of  the  symmetry  about  the 
z axis  of  the  contact  surfaces  in  figure  5.1,  the  pressure  distribution 
must  also  be  symmetric  about  z.  Hence,  p (x,  3)  may  be  replaced  by 
p (x)  in  equation  (5.7),  (see  fig.  5.3)  and  p ($,  3)  in  equation 
(5.8)  may  be  replaced  by  p (<}>).  Making  these  substitutions  and 
combining  equations  (5.3),  (5.7)  and  (5.8)  the  contact  criterion, 
equation  (5.3)  becomes 


S - f £44  f ~ S + 


♦ /xx  A /?, 2 / t(Y)G, (XPX  y'dpJY' 

JJL, 


(5.9) 


The  boundary  conditions  require  that  the  separation,  S , be  zero  with 
positive  pressure  inside  the  contact  regions  ^ and  fig.  Also,  the 
separation  must  be  positive  with  zero  pressure  outside  the  contact 
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Fig.  5.3.  Axisymmetric  pressure  loading  on  sphere 
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regions.  In  symbolic  terms 


ltJ$U>E  U1 


(5.10a) 


?( 4>),  t(Y)-0  iNsioe  j\ 


(5.10b) 


Sv  >0 


Outside  si 


(5.10c) 


fC^),  f(Y)  - O Oursme  n 


(5. lOd) 


It  is  required  to  determine  the  contact  regions,  ^ and  the 

pressure  distributions  p (ip ) and  p ( <t> ) , and  the  approach,  6,  such 

that  relationships  (5.9)  and  (5.10)  are  satisfied. 

The  determination  of  the  contact  regions,  ^ and  poses  a 

major  problem  in  the  solution  of  equation  (5.9).  However,  choosing 

★ * 

some  tentative  "candidate"  regions,  ft-j  and  ^ * will  establish 
equation  (5.9)  as  an  integral  equation  of  the  first  kind  which  can 
be  solved  using  the  "Simply  Discretized"  method  of  Singh  and  Paul 
[1974],  In  the  conformal  contact  problem  of  a sphere  indenting  a 
spherical  seat,  it  is  known  a priori,  that  the  boundary  of  ft,  will 
always  be  a circle  on  the  sphere  defined  by  . Also,  the 

Ula  X 


boundary,  <b  , of  ft9  may  be  chosen  such  that  s9  = s,  for  the 
max  l.  l. 


* 


* 


boundaries  of  ft-|  and  • Thus 


(5.11) 
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The  region  ^ is  written  without  an  asterisk  since  it  can  be  defined 

as  the  exact  contact  region  on  the  sphere  corresponding  to  some 

* 

unknown  force  F.  ^ 1S  denoted  with  the  asterisk  since  it  is  chosen 

according  to  equation  (5.11).  Further  details  on  the  refinement  of 

* 

the  tentative  region  ^ are  discussed  later. 

A "Simply  Discretized"  solution  is  found  by  subdividing  ft-j 

* 

and  ^ into  a large  number  of  small  cells.  The  normal  pressure  is 

replaced  by  a piecewise  constant  pressure  field  (pressure  p^  and 

* 

cell  i).  Thus  and  are  divided  into  N cells  apiece,  such  that 
the  centroid  of  cell  i of  fi-j  merges  with  the  centroid  of  cell  i of 

■k 

^2  • This  is  achieved  by  first  choosing  the  cells  on  fi-j  and  then 

using  the  relation  S2  = s.j  to  determine  the  corresponding  cells  on 
* 

!w!  ^ • 

Because  of  the  known  symmetry  of  p (c|>)  and  p (ip) , the  cells 

on  each  surface  will  be  chosen  as  rings  symmetric  about  the  n-|  axis 

as  shown  in  figure  5.4  for  the  sphere.  Cell  i will  be  located  on  the 

sphere  between^  and  i|k+.j  while  that  on  the  seat  will  be  between  <£.. 

and  4>.j+.| . N field  points  are  chosen  on  each  surface  such  that  the 

"t  h t h 

location  of  the  j field  point  is  within  the  j n cell.  Equation 
(5.9)  may  be  written  in  discretized  form  for  field  point  j as, 

$h  - f ^c&4  t*  “ ^ cC  t 
+«*xzpifG,(ri,p),  y;  f\ d„  e^r^'dr 

^ li 

Z Pi  \ Gi  (tj,  l) t,  'A  fid f (5.12) 
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A 


Fig.  5.4.  Axi symmetric  cells  on  sphere 
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The  integrals  in  equation  (5.12)  may  be  evaluated  numerically.  The 

location  of  all  the  field  points  will  be  along  the  contour  curve 

located  in  the  n-|  - Y-|  and  P^anes  so  that  B'j  will  always  be 

considered  zero.  Equation  (5.12)  may  be  written  for  each  field  point. 

Thus  N linear  algebraic  equations  are  generated  in  N + 1 unknown 

(N  values  of  P^,  and  6).  One  additional  equation  is  needed  to  produce 

a unique  solution  of  the  unknown  variables.  The  last  equation  is 

generated  by  writing  equation  (5.12)  for  one  additional  field  point. 

The  location  of  this  field  point  is,  in  theory,  arbitrary  but,  as  will 

be  discussed  later,  the  results  are  sensitive  to  this  choice. 

Having  generated  N + 1 linear  algebraic  equations  in  N + 1 

unknowns  the  piecewise  constant  pressure  distribution,  P. , and  6 may 

be  found.  It  now  remains  to  check  the  validity  of  the  method  used  to 
★ 

determine  and  the  location  of  two  points  A and  B which  merge  after 

deformation.  In  each  case,  the  merging  points  or  boundaries  were 

initially  chosen  such  that  s^  = Sy  It  will  now  be  illustrated  how  the 

"point-mating"  procedure,  described  in  chapter  3 may  be  utilized  to 

★ 

refine  the  choice  of  the  outer  boundary  of  and  of  points  on  the 
seat  such  that  they  merge  with  the  appropriate  points  on  the  sphere. 

The  total  separation  after  deformation  of  points  A and  B may 
be  computed  from  equation  (3.13a).  For  the  case  of  a sphere  indenting 
a spherical  seat,  the  value  of  SQ  is  zero.  Furthermore,  the  contact 

fc) 

criterion  required  that  is  zero.  Therefore,  the  only  non-zero 
component  of  separation  is  St  where 


a 604.  a (ot,- ux)  + ft  - 


(5.13a) 
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and 


64-6  trt  (tt-cx) 

ft  r / f I M&\  f 


(5.13b) 


(5.13c) 


The  components  of  displacements  u-j  and  may  be  computed  using  the 
pressure  distribution  (P^)  of  the  "Simply  Discretized"  solution  by 


where  and  are  the  influence  functions  for  displacements  in  the 
tangential  direction  on  the  sphere  and  seat  respectively.  is 
defined  by  H in  equation  (4.24)  and  is  defined  by  H in  equation 
(4.32).  Having  evaluated  S^.,  at  all  field  points,  the  "point-mating 
procedure"  may  be  used  to  find  coordinates  for  a new  set  of  field 
points  on  the  seat  which  upon  repeating  the  solution  procedure  just 
described  will  give  rise  to  new  values  of  St  which  are  smaller  than 
those  previously  calculated.  Originally,  mating  points  on  the  sphere 
and  seat  were  located  by  assuming  that  the  distances  along  their 
respective  contour  curves  were  equal , i .e. , s^  = s-j  as  illustrated  in 


(5.14a) 


and 


)Autjb‘ df'd 


(5.14b) 


I 
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figure  3.1.  The  "point-mating  procedure"  discards  this  relation  after 
the  initial  solution  and  replaces  it  with  relations  (3.16  a,b). 

Repeated  iterations  yield  values  of  which  become  smaller  if  the 
process  converges.  The  limits  on  the  final  values  of  depend  largely 
on  the  cell  density  and  the  values  of  Numerical  experiments 

have  shown  that  for  large  values  of  ^max  (>  10°),  with  few  cells 
(<  10),  the  final  value  of  which  can  be  achieved  is  of  the  order 

of  e (w-j  + W£)  where  e = 0.1.  For  problems  with  up  to  15  cells  values 
of  e = 0.01  can  be  achieved.  Considering,  essentially  Hertzian 
problems  (^max  < .1°),  e = 0.01  may  be  achieved  on  the  first  solution. 
The  location  of  (f>max  defining  may  a^so  be  refined  using  the 
"point-mating  procedure."  In  most  cases  this  was  performed  auto- 
matically since  4>max  was  chosen  as  the  additional  field  point. 

Having  generated  the  displacements  w-j , w^,  u-j  and  u^  via 
equations  (5.7),  (5.8)  and  (5.14)  at  field  points  along  a contour 
curve,  the  strains  cOD,  e. , on  the  sphere  and  cOD  and  e. . on  the  seat 

PP  w PP  <p<p 

may  be  formulated.  These  quantities  will  be  needed  later  for 
accuracy  analysis.  It  can  be  shown^  that  the  strains  on  the  surface 
of  a sphere  for  the  axi symmetric  set  of  displacements  are 


+ 


u,  (*) 

R, 


Mm 


(5.15) 


^okolnikoff,  2d  ed.,  1956,  p.  184. 
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and 


fn>.  \ . ! dUi('V'i)  cOiO'i) 

l)~  R,  dy  R,  (5' 

J_  I- 

where  i|;.  is  the  coordinate  of  the  l field  point.  Similarly,  for 
the  spherical  seat 


Rl 


tot  (${) 


(5.17) 


and 


/ dUi(fr) 
Rt  ol? 


4 


Wt  (fii) 

Ri 


(5.18) 


±.  L 

where  (|>.  is  the  coordinate  of  the  iL  field  point  on  the  seat.  These 
strains  may  be  computed  using  finite  difference  approximations  for  the 
derivatives  in  equations  (5.16)  and  ( 5 . 1 8 ) 1 . 

5.3  Numerical  Procedures 

In  each  "Simply  Discretized"  solution  the  boundaries  on 
and  are  defined  by  to  and  <j)  respectively.  Each  contact  region 

c.  MTaX  max 

was  divided  into  N cells  by  first  subdividing  ib  and  <t>  into  N 

3 max  rmax 

equal  intervals  subtending  equal  arc  lengths.  The  cells  were  then 
^Carnahan,  [1969]  p.  431  „ 
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defined  by  the  surface  generated  by  revolving  the  arc  lengths  about 
the  n-j  and  n2  axis.  Typical  cells  on  the  sphere  are  illustrated  in 
figure  5.4.  All  field  points  were  located  on  the  contour  lines 
defined  by  the  n-|  - y-|  and  n2  - Y2  planes,  (see  fig.  5.4) 

The  integrals  in  equations  (5.12)  and  (5.14a,  b)  were 
evaluated  in  part  by  analytic  means.  Consider  the  integration  of  an 
influence  function  corresponding  to  the  i^  field  points  at  ^ (<J>^) 
over  the  region  (^2j.).  For  all  cases  where  i f j,  the  integra- 
tion was  performed  by  Gaussian  quadrature.  When  i = j,  i.e.,  when 
the  field  point  is  located  within  the  cell  of  integration,  the 
integrand  is  singular  within  the  region  of  integration  and  the 
singularity  is  located  at  the  field  point.  Integration  over  the 
singularity  is  performed  analytically  while  that  over  the  remaining 
portion  of  the  cell  is  performed  numerically  using  Gaussian 
quadrature. 

Consider  the  portion  of  a cell  on  a sphere  of  radius  R-j  near 
the  field  point  at  (^  , 0)  as  shown  in  figure  5.5.  The  region  is 
bounded  by  arcs  defined  by  3 = ±A,  tp  = ipc  - A and  tp  = ij>  + A 
(A  < < 1).  Similarly,  consider  a region  on  the  seat.  For  both  the 
sphere  and  seat  when  small  cells  are  used,  the  boundaries  of  the  cell 
fall  within  the  limits  of  this  region.  In  this  case  the  values  of 
i|k  and  (}).  defining  the  cell  boundaries  were  used  to  bound  the  region. 
For  small  A,  taken  in  this  analysis  to  be  1/4  degree,  the  region 
around  the  singularity  on  the  sphere  will  be  approximated  by  a small 
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Fig.  5.5.  Small  region  surrounding  a field  point  on  a 
spherica.1  surface. 


Fig.  5.6.  Planar  approximation  to  a region  surrounding  a 
field  point  on  a spherical  surface. 


^ > 


I 


no 

planar  annular  element  with  inner  radius  r^ , outer  radius  r^,  and 
half  angle  p where 

| 

A}  - R,  Aw  ( yc  - Yc  (5.19a) 

>1  - R\  a™  (Yc*A)  / <004.  Yc  (5.19b) 

i 

and 

I 

J*  = A /C^i  Tc  (5.19c) 

Similar  relations  in  terms  of  cj>  , and  A may  be  written  for  the 
definition  of  r-j , r^  and  p for  the  seat.  The  derivation  of  r^ , r^ 
and  p is  shown  in  appendix  P.  The  displacement  at  the  field  point 
C due  to  a constant  force  applied  to  a general  point  A within  the 
region  is  governed  by  the  Boussinesq  influence  function  as  shown  in 
appendix  A.  The  integration  of  the  Boussinesq  influence  function 
over  an  annular  element  is  derived  in  appendix  G.  Therefore,  the 
integration  over  a cell  i for  field  point  i was  performed  in  two 
parts.  First  a small  annular  element  surrounding  the  singularity 
(field  point)  was  defined,  and  the  integral  within  that  region  was 
computed  using  the  analytical  solutions  in  appendix  G.  Second  the 
remaining  portion  of  the  integral  was  computed  using  2-dimensional 
10  point  Gaussian  quadrature.  Similar  approximations  were  made  for 
the  spherical  seat. 

The  "point-mating  procedure"  was  employed  to  identify  the  set 
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of  field  points  on  the  seat  which  merge  with  the  field  points  on  the 
sphere.  The  iteration  was  terminated  when  values  of  S for  all  field 
points  were  within  their  respective  values  of  [ (w-j  + w^)  e ] ^ . e was 
taken  to  be  0.1  for  preliminary  results  while  accurate  solutions  were 
obtained  by  choosing  e to  be  0.01. 

For  all  cases  that  follow,  the  extra  field  point  was  located 

at  the  outer  boundary  of  the  contact  region.  With  this  location,  the 

final  determination  of  <f>  on  the  seat  is  automatically  performed, 

ma  x 

i.e.,  4>max  is  the  final  coordinate  of  the  extra  field  point.  The 
"Simply  Discretized"  method  of  solution  always  provided  stable 
pressure  distributions  even  though  Singh  and  Paul  [1973]  noted  that 
for  the  problems  they  treated,  this  method  generated  an  ill- 
conditioned  set  of  equations  which  resulted  in  unstable  (widely 
varying  and  negative)  pressure  distributions.  This  stability  can 
possibly  be  attributed  to  the  choice  of  axi symmetric  cell  distribu- 
tions which  had  not  previously  been  attempted  in  the  work  of  Singh 
and  Paul.  The  solutions  will  be  termed  "quasi  stable"  because  the 
solution  becomes  unstable  if  the  location  of  the  extra  field  point  is 
moved  within  the  contact  region.  The  success  achieved  when  the  extra 
field  point  is  located  on  the  boundary  of  the  contact  region  can 
possibly  be  attributed  to  the  fact  that  the  location  is  farthest  from 
all  other  field  points,  on  the  average,  than  any  location  inside  the 
contact  region.  This  choice  could  yield  the  most  independent  extra 
equation  thus  providing  a more  stable  solution. 

Having  generated  N + 1 equations  in  N + 1 unknowns 
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(N  values  of  and  6)  as  described  in  section  5.2,  the  set  was  first 

reduced  to  N equations  in  N unknowns  by  subtracting  the  equation 

written  for  the  extra  field  point  from  the  others  thus  eliminating  6 

from  the  set.  The  remaining  N equations  were  solved  using  Gaussian 

s t 

elimination.  The  approach,  6,  was  then  computed  from  the  N + 1 
equation. 

The  strains  corresponding  to  equations  (5.17),  (5.18),  (5.19) 
and  (5.20)  were  computed,  from  equations  (5.15)  - (5.18),  using  the 
values  of  w^ , w^,  u^  and  u^  at  each  field  point  on  the  sphere  and  seat. 
The  required  derivatives  du-j/dij;  and  du^/dtj)  were  evaluated  using 
central  difference  formulae  for  field  points  i where  i / 1 or  N.  For 
i = 1 forward  difference  formulae  were  used  while  for  i = N backward 
difference  formulae  were  usedJ 

The  total  force  applied  to  the  sphere  can  be  calculated  from 
the  pressure  distribution  found  in  the  analysis.  Consider  the  discre- 
tized pressure  distribution  P..  on  cell  i.  The  force  in  the  - n-j 

2 

direction  from  P.  applied  to  a small  sector  area  2tt  sin  ipdip  is 
2 

2ttR^  P^  cos\p  sin  ipdip.  Integrating  this  between  ^ andi|A+^,  the  force 

2 2 

in  the  - n-j  direction  due  to  P^  on  cell  i is  ttR^  P.  • (sin  ip^+-|  - 
2 

sin  ip.).  Thus  the  total  force  applied  to  the  sphere  can  be  computed 
as 

F - Z irfi*  % (5  zo) 

4*3/ 

^See  Carnahan  [1969]  p.  431. 
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5.4  Numerical  Results 

A computer  program  termed  CONSPHERE  was  written  to  analyze 
the  conformal  contact  between  an  elastic  sphere  and  seat.  The 
following  numerical  example  was  considered 

R-|  =1  in 

R^  = 1.01  in 

V1  = v2  = ^-25 

E]  = E2  = 30  x io6  psi 

'I'max  = °'5  de9r^es  (5.21) 

In  this  example  and  the  ones  to  follow,  the  value  of  Poisson's 
ratio  was  chosen  to  be  0.25  so  that  the  results  could  be  compared  to 
those  of  Goodman  and  Kerr  [1965].  The  results  are  presented  in 
dimensionless  form.  Let, 


<?>?,  Rx 

R x ~ R i 


(5.22) 


and 


(5.23) 


114 


Then  define 


Dimensionless  Pressure  in  cell  i,  P.  = P.k 

i i 


(5.24a) 


* R 

Dimensionless  Load,  F = F ^“2~ 

m 


(5.24b) 


Ri 

* I 


Dimensionless  Radii,  R,  = r 

i 


m 


(5.24c) 


* Ro 

"2 

m 

* 5 

Dimensionless  Approach,  6 = p 


m 


(5.24d) 


(5.24e) 


Dimensionless  Ratio  of  Radii  = 


R, 

r: 


(5. 24f ) 


The  contact  region  was  divided  into  15  cells  such  that 
^j+1  = 0.033  degrees.  The  pressure  distribution  obtained  from 

CONSPHERE  is  compared  in  figure  5.7  to  the  pressure  distribution 

predicted  by  Hertzian  theory  for  the  same  applied  force  of 

* -14  * -8 

F = 0.7294  x 10  . The  approach  6 was  found  to  be  0.3858  x 10  in, 

★ -8 

while  the  value  predicted  by  Hertzian  theory  is  6 = 0.3895  x 10  in. 

This  problem  falls  within  the  domain  of  Hertzian  theory,  and  the 
comparison  of  these  results  indicates  that  the  solution  produced  by 
CONSPHERE  is  in  general  agreement  with  those  of  Hertz. 

Now  consider  the  non-Hertzian  problem  where  ^max  = 30  degrees 
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and  the  dimensions  and  elastic  constants  are  the  same  values  as  used 

l 

in  the  previous  example.  Figure  5.8  illustrates  the  pressure  distri- 
bution between  the  sphere  and  seat.  Also  shown  is  the  Hertzian 
pressure  distribution  for  the  same  applied  force.  It  should  be  noted 
that  the  conformal  contact  solution  predicts  a higher  peak  stress  and 

a smaller  contact  region  than  the  Hertzian  solution.  The  approach 

* -4 

for  the  conformal  solution  is  6 = 0.1281  x 10  and  the  total 

* -7 

compressive  force  was  found  to  be  F = 0.1780  x 10  . The  radial  and 

tangential  displacements  are  tabulated  in  tables  5.1  and  5.2.  It  is 
interesting  to  note  that  the  tangential  displacements  on  the  sphere 
are  all  positive  (increasing  \p  is  positive  direction)  while  those  on 
the  seat  are  negative  (increasing  <p  is  positive).  This  can  be  under- 
stood if  one  considers  the  sphere  to  be  flattened  out  while  the  seat 
is  a depression  which  elongates  or  grows  deeper. 

Additional  problems  were  solved  in  order  to  compare  the  load- 
approach  relationship  to  that  obtained  with  the  solutions  of 
Hertz  [1881]  and  Goodman  and  Keer  [1965].  Figure  5.9  illustrates  the 
load-approach  curves  for  Hertzian  theory,  the  theory  of  Goodman  and 
Keer  and  experimental  data  reported  by  Goodman  and  Keer.  The  results 
are  plotted  for  half  angles  of  contact  between  0°  and  20°.  Figure  5.9 
clearly  illustrates  a strong  correspondence  between  the  present 
theory  and  the  experimental  data  reported  by  Goodman  and  Keer.  The 
load-approach  curve  for  the  Hertzian  theory  indicates  more  compliance 
than  that  of  the  other  theories  while  the  load-approach  relation  of 
the  Goodman  and  Keer  theory  is  less  compliant  than  the  others.  All 


SOLUTION  BY  CONSPHPRS 


117 


O 

£ 


ki 

u 

$ 

V) 

5 


o 

O 

ro 


X 

fO 


Fig.  5.8.  Pressure  distribution  between  sphere  and  seat. 
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TABLE  5.1 

SURFACE  DISPLACEMENTS  ON  SPHERE 


Ip 

[Deg] 

Radial 

Displacements 
w-j  x lo4  in 

Tangential 

Displacements 

u-j  x 104  in 

1.0 

14.12 

0.03869 

3.0 

14.04 

0.1194 

5.0 

13.87 

0.2035 

7.0 

13.62 

0.2841 

9.0 

13.29 

0.3578 

11.0 

12.88 

0.4300 

13.0 

12.39 

0.4968 

15.0 

11.82 

0.5580 

17.0 

11.16 

0.6123 

19.0 

10.43 

0.6635 

21.0 

9.627 

0.7116 

23.0 

8.745 

0.7562 

25.0 

7.792 

0.7998 

27.0 

6.766 

0.8475 

29.0 

5.672 

0.9165 

30.0 

5.099 

0.9711 

iji  = 30°,  R,  = 1.00  in,  R0  = 1.01  in, 
m x I c 

E^  = E^  = 30  x 10^  psi , vi  = v2  = 0-25 
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TABLE  5.2 

SURFACE  DISPLACEMENTS  ON  SPHERICAL  SEAT 


[Dig] 

Radial 

Di splacements 
w 2 x 104  in 

Tangential 

Displacements 

x 104  in 

0.9888 

11.73 

- 0.1789 

2.966 

11.66 

- 0.5307 

4.944 

11.53 

- 0.8735 

6.922 

11.32 

- 1.211 

8.899 

11.04 

- 1.542 

10.88 

10.70 

- 1.858 

12.85 

10.29 

- 2.158 

14.83 

9.812 

- 2.439 

16.81 

9.267 

- 2.699 

18.79 

8.654 

- 2.929 

20.77 

7.974 

- 3.125 

22.74 

7.228 

- 3.286 

24.72 

6.414 

- 3.404 

26.70 

5.534 

- 3.470 

28.68 

4.587 

- 3.463 

29.67 

4.089 

- 3.418 

^max  = 30°’  R1  = 1,00  in’  R2  = 1,01  in’ 
E-|  = E^  = 30  x 10^  psi,  v-j  = v2  = 0-25 
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□ - CONSPHERE  (v=0 .25 ) 

O - CONSPHERE  (v=0.30) 

A - EXPERIMENTAL  (v=0.33) 


Fig.  5.9  Load-approach  relationship  for  conformal  contact  of 
a sphere  and  spherical  seat. 
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2 

of  the  theories  agree  for  loads  less  than  F (1  - v)  / R-j  E = 2 while 
the  experimental  data  near  this  region  deviates  from  the  theories. 
Goodman  and  Keer  [1965]  attribute  this  discrepancy  to  experimental 
error. 

f 

The  effect  of  Poisson's  ratio  on  the  load-approach  curve  was 
also  studied.  The  variation  of  force  and  approach  for  values  of 
Poisson's  ratio  of  v = 0.25  and  v = 0.30  were  plotted  together  in 
figure  5.9.  With  increasing  Poisson's  ratio  the  materials  were  found 
to  be  less  compliant. 

A comparison  was  made  of  the  val ues , predicted  by  the  various 
theories,  of  the  radius  of  the  contact  region  a = sin  ( and  6 . 
Figure  5.10  illustrates  the  results  of  Hertz,  Goodman  and  Keer,  and 
the  present  theory.  The  results  from  CONSPHERE  fall  much  closer  to 
the  Hertzian  theory  than  those  of  Goodman  and  Keer. 

Knowledge  of  the  displacements  at  discrete  points,  namely, 
the  field  points,  enables  one  to  also  calculate  the  surface  strains 
given  in  equations  (5.15)  - (5.18).  Both  the  displacements  and 
strains  for  several  problems  will  be  compared  to  those  obtained  in 
the  analysis  of  Goodman  and  Keer  [1965].  It  will  be  shown  through 
this  comparison  that  the  assumptions  used  by  Goodman  and  Keer  can 
produce  erroneous  displacements  and  strains  in  problems  where  the 
contact  angle  exceeds  i b = 60°. 

In  the  derivation  of  the  contact  criterion  used  by  Goodman 
and  Keer  (see  appendix  L)  it  is  assumed  that  points  on  both  bodies, 
which  are  initially  equidistant  from  the  axis  of  symmetry  come  into 
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conformal  contact  of  sphere  and  spherical  seat. 
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contact  after  deformation.  To  show  that  this  assumption  can  not  be 

true  for  large  angles  of  contact,  consider  the  contact  of  an  elastic 
sphere  in  a rigid  seat,  as  shown  in  figure  5.11.  According  to 
Goodman  and  Keer  points  A and  B which  merge  after  deformation,  are 
located  such  that  R-j  sin  ip  = R?  sin  <J>.  Since  the  seat  is  rigid, 
point  A on  the  sphere  will  merge  with  point  B on  the  seat  solely  due 
to  a displacement  uz  of  point  A.  This  displacement  may  be  viewed  as 
having  components  in  the  r and  t directions,  i.e.,  w-j  and  u^ , where 
r is  directed  radially  inward  on  the  sphere  and  t is  perpendicular 
to  r,  i . ,e. , 


fU,  t 


(5.25) 


where  (see  appendix  M for  derivation): 


l(Tt/  s Rt  (l-AO*#) + 6 (5.26) 

IV,  - / Wj  / £04,  'Y  (5.27) 


Wz  - ItitlAiH  'f 


(5.28) 
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Fig.  5.11.  Sphere  in  contact  with  seat  (displacements  of 
points  A and  B constrained). 
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and 

- CM4MTI  ( ^ Aon  ^ ) (5.29) 

Now  consider  the  case  where 

R.|  = 1.0  in. 

R2  = 1.001  in. 

E]  = 30  x lo6  psi,  E2  = 30  x io10  psi 

v-|  = 0.25,  v2  = 0.25 

^max  ‘ 60°  <5-30> 

Since  E2  >>  E-j  the  seat  will  be  considered  completely  rigid 
and  equations  (5.27)  and  (5.28)  will  be  used  to  compute  w^  and  u-j 
for  results  obtained  in  the  Goodman  and  Keer  analysis.  The  displace- 
ments for  the  present  conformal  theory  of  section  5.2  were  computed 
using  equations  (5.7),  (5.8)  and  (5.14).  The  displacement  results 
for  the  problem  at  hand  are  compared  in  tables  5.3  and  5.4.  Also 
tabulated  are  the  displacements  that  are  obtained  when  the  pressure 
distribution  obtained  through  CONSPHERE  was  applied  to  a finite 
element  model  of  the  sphere.^  The  radial  displacements  of  all 
solutions  agree  well.  The  tangential  displacements  obtained  through 
the  Goodman  and  Keer  analysis  are  much  higher  than  those  of  CONSPHERE 

Vhe  paper  of  Goodman  and  Keer  [1965]  did  not  provide  any  data 
for  angles  of  contact  above  20  degrees.  Therefore  no  comparable 
finite  element  model  could  be  analyzed  using  their  pressure  distri- 
bution. 
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TABLE  5.3 

COMPARISON  OF  RADIAL  DISPLACEMENTS  w] 


[Deg] 

Goodman 
and  Keer 

w-j  x 103  in 

CONSPHERE 
w-j  x 103  in 

Finite 

Element 

w-j  x 103  in 

3.0 

1.171 

1 .166 

1.184 

9.0 

1.147 

1.142 

1.157 

15.0 

1.100 

1.095 

1 .109 

21 .0 

1.029 

1.025 

1.040 

27.0 

0.9368 

0.9330 

0.9480 

33.0 

0.8231 

0.8197 

0.8315 

39.0 

0.6895 

0.6864 

0.6965 

45.0 

0.5374 

0.5348 

0.5364 

51.0 

0.3686 

0.3663 

0.3535 

57.0 

0.1850 

0.1828 

0.1738 

= 60°,  R,  = 1.000  in,  R0  = 1.001  in, 
rmax  1 2 

E]  = 30  x IQ6  psi,  E2  = 30  x to10  psi , v = v2  = 0.25 


127 


TABLE  5.4 

COMPARISON  OF  TANGENTIAL  DISPLACEMENTS  u] 


4> 

[Deg] 

Goodman 
and  Keer 
ln4  . 

u-j  x 10  in 

CONSPHERE 
,n4  • 

u-j  x 10  in 

Fi  ni  te 
Element 

4 

u-j  x 10  in 

3.0 

0.6135 

0.2001 

0.2065 

9.0 

1.816 

0.5972 

0.6160 

15.0 

2.946 

0.9737 

1 .002 

21.0 

3.951 

1 .306 

1.338 

27.0 

4.773 

1.576 

1 .607 

33.0 

5.345 

1.770 

1.806 

39.0 

5.584 

1 .881 

1.910 

45.0 

5.374 

1 .892 

1 .915 

51.0 

4.552 

1.811 

1 .858 

57.0 

2.849 

1 .657 

1 .718 

^max  = 60°’  R1  = in,  = l.OOl  in, 

E1  = 30  x 106  psi,  E2  = 30  x 1010  psi , v]  = v2  = 0.25 


128 


or  the  finite  element  solution. 

Having  computed  the  displacement  fields,  the  strains  eOD  and 

38 

may  be  calculated  using  relations  (5.15  - 5.18).  They  are 
tabulated  in  tables  5.5  and  5.6  for  the  analysis  by  Goodman  and  Keer, 
CONSPHERE  and  finite  element  theory.  It  can  be  seen  that  the 

88 

strains  are  zero  for  the  analysis  of  Goodman  and  Keer  which  is 
expected  since  the  circle  defined  by  ip  does  not  enlarge  or  shrink 
after  deformation.  The  strains  e from  CONSPHERE  are  in  general 
agreement  with  those  produced  through  finite  element  analysis.  The 

computation  of  e reveals  that  near  the  boundary  of  the  contact 

W 

region,  those  strains  predicted  by  the  Goodman  and  Keer  analysis  are 
much  larger  than  the  solutions  of  CONSPHERE  arid  the  finite  element 
analysis.  Furthermore,  the  results  of  CONSPHERE  agree  with  those 
of  the  finite  element  analysis. 

Both  p and  e,  . decrease  as  \p  increases  according  to 

88  # 

CONSPHERE  and  the  finite  element  analysis.  This  is  to  be  expected 
since  the  pressure  diminishes  when  ip  increases.  On  the  other  hand 

the  values  of  e,  , predicted  by  the  Goodman  and  Keer  model  increase  as 

W 

ip  increases.  Finally,  the  values  of  enn  and  e,.  are  nearly  equal  for 

38  # 

small  ip  in  the  s'olution  of  CONSPHERE  and  finite  element  analysts 
which  is  expected  in  this  axi symmetric  case. 

5.5  Conclusions 

The  problem  of  a sphere  indenting  a spherical  seat  has  been 
solved.  The  pressure  distribution  for  ip  = 0.50°  has  been  shown  to 
compare  closely  with  the  Hertzian  solution,  (fig.  5.7)  For  to  = 30°, 

MlaX 
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TABLE  5.5 

CONPARISON  OF  STRAINS  e. . 


Ip 

[Deg] 

Goodman 
and  Keer 

x 103  in/ in 

CONSPHERE 
x 103  in/ in 

Fini te 
El  ement 

x 103  in/ in 

3.0 

- 0.002744 

- 0.7868 

- 0.7894 

9.0 

- 0.02506 

- 0.7730 

- 0.7692 

15.0 

- 0.07172 

- 0.7567 

- 0.7563 

21.0 

- 0.1472 

- 0.7377 

- 0.7355 

27.0 

- 0.2593 

- 0.7115 

- 0.7121 

33.0 

- 0.4210 

- 0.6741 

- 0.6626 

39.0 

- 0.6545 

- 0.6280 

- 0.6181 

45.0 

- 0.9975 

- 0.5681 

- 0.5437 

51.0 

- 1.520 

- 0.4788 

- 0.4171 

57.0 

- 2.360 

- 0.3299 

- 0.2873 

iL  = 60°,  R,  = 1.000  in,  R0  = 1.001  in, 

FTla  X I c 

E]  = 30  x 106  psi,  E2  = 30  x 1010  psi,  v]  = v2  = 0.25 
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TABLE  5.6 

COMPARISON  OF  STRAINS  eoc 

pp 


[Deg] 

Goodman 
and  Keer 

in/ in 

CONSPHERE 
x 103  in/ in 

Finite 

Element 

x 103  in/ in 

3.0 

0.0 

- 0.7842 

- 0.7900 

9.0 

0.0 

- 0.7653 

- 0.7680 

15.0 

0.0 

- 0.7318 

- 0.7354 

21.0 

0.0 

- 0.6849 

- 0.6912 

27.0 

0.0 

- 0.6237 

- 0.6325 

33.0 

0.0 

- 0.5471 

- 0.5535 

39.0 

0.0 

- 0.4542 

- 0.4604 

45.0 

0.0 

- 0.3455 

- 0.3444 

51.0 

0.0 

- 0.2196 

- 0.2029 

57.0 

0.0 

- 0.07525 

- 0.06036 

^max  = 60°’  R1  = K00°  in’  R2  = 1-001  in’ 

E]  = 30  x io6  psi,  E2  = 30  x 1010  psi,  v]  = = 0.25 
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the  resulting  pressure  distribution  has  a higher  peak  stress  and 
smaller  contact  area  than  the  Hertzian  solution  for  the  same  applied 
force.  The  load-approach  curve  was  plotted  and  was  found  to  agree 
closely  with  the  experimental  results  reported  by  Goodman  and  Keer. 

The  plot  of  approach  vs.  contact  area  indicates  that  the  results  are 
nearly  equal  to  those  of  the  Hertzian  theory  for  angles  up  to  20°. 

The  displacement  field  and  resulting  strains  were  found  to  be 
reproduced  when  a finite  element  analysis  was  made  of  the  sphere  under 
the  pressure  field  predicted  by  CONSPHERE  for  \p  = 60°  (the  seat 
being  rigid).  It  was  shown  that  in  using  the  constraint  imposed  by 
Goodman  and  Keer  in  their  analysi s,  larger  tangential  displacements 
and  strains  will  be  produced.  In  view  of  the  fact  that  the 
uniqueness  theorem  for  contact  theory  (proven  by  J.  J.  Kalker  1971) 
guarantees  a unique  pressure  field  and  displacement  field  for  a given 
contact  area,  it  can  be  concluded  that  the  total  solution  predicted 
by  CONSPHERE  is  correct  since  the  finite  element  model  reproduces  the 
same  displacement  field  as  CONSPHERE  when  subjected  to  the  interfacial 
pressure  predicted  by  CONSPHERE. 


6.  CONFORMAL  ELASTIC  CONTACT  OF  A CYLINDER 
INDENTING  A CYLINDRICAL  CAVITY 

6.1.  Introduction 

The  solution  of  the  two  dimensional  contact  problem  of  two 
cylinders  in  contact  or  the  problem  of  a cylinder  indenting  a 
cylindrical  seat  can  be  obtained  from  Hertzian  theory  by  allowing  the 
radii  of  curvature  of  each  body  to  become  infinite  in  one  direction.^ 
However,  such  a solution  is  only  valid  within  the  assumptions  of 
Hertz,  i.e.,  the  in  plane  dimensions  of  the  contact  area  must  remain 
small  compared  to  the  in  plane  radii  of  curvature.  Therefore,  Hertz's 
solution  is  not  appropriate  for  moderate  loads,  when  the  difference  in 
the  radii  of  the  cylinder  and  cylindrical  seat  is  small. 

A more  recent  theory,  pertaining  specifically  to  the  problem 
of  a cylinder  indenting  a cylindrical  seat  has  been  published  by 
Sjtaerman  [1949]  for  the  problem  where  the  contact  pressure  on  the 
cylinder  and  seat  are  equilibrated  by  identical  pressures  located  at  tt 
radians  from  the  contact  region.  Sjtaerman's  solution  is  based  on  the 
formulation  of  a contact  criterion  in  the  radial  direction  of  a polar 
coordinate  system  fixed  at  the  center  of  the  cylinder.  He  does  not 
consider  displacements  tangential  to  the  surface.  In  order  to  compute 
the  radial  displacements  within  the  contact  area,  Sjtaerman  forms  the 

^Timoshenko  and  Goodier,  3d.  ed.  [1970],  pp.  418-20. 
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integral  of  the  unknown  pressure  distribution  times  the  influence 
function  for  the  cylinder  and  seat.  These  integrals  are  incorporated 
into  the  contact  criterion  and  the  unknown  pressure  field  is  determined 
using  a finite  difference  technique. 

A third  theory  specific  to  the  problem  of  a disc  contacting  a 
hole  in  an  infinite  plate  has  been  published  by  Persson  [1964].  In 
contrast  to  the  work  of  Sjtaerman  and  the  present  analysis,  Persson 
considers  the  disc  to  be  equilibrated  by  a force  located  at  the  center 
of  the  disc  with  the  seat  being  fixed  at  infinity.  Persson  initially 
assumes  the  existence  of  both  tangential  and  radial  displacements  of 
surface  points.  He  found  that  by  assuming  the  contact  region  to  be 
circular  and  neglecting  second  order  quantities,  the  contact  criterion 
was  independent  of  tangential  displacements  and  is  identical  to  that 
derived  by  Sjtaerman.  Persson  proceeds  to  develop  the  final  form  of 
the  contact  criterion  in  terms  of  a singular  integro-differential 
equation  which  he  solves. 

The  problem  of  an  elastic  cylinder  indenting  an  elastic  cy- 
lindrical seat  is  solved  in  this  chapter  using  the  conformal  theory 
developed  in  chapter  3.  No  assumptions  pertaining  to  the  tangential 
displacements  are  retained  in  the  final  solution.  The  loading  condi- 
tions applied  are  the  same  as  those  used  in  the  Sjtaerman  analysis. 

Section  6.2  contains  the  formulation  of  the  problem  and  the 
numerical  procedures  are  discussed  in  section  6.3.  The  results  are 
compared  to  those  of  Hertz,  Sjtaerman  and  Persson  in  section  6.4. 
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6.2  Formulation 

Consider  the  conformal  contact  of  an  elastic  cylinder  of 
radius  R-j  indenting  an  elastic  cylindrical  seat  of  radius 
(R^  > R-j).  It  is  assumed  that  the  bodies  are  equilibrated  by  pressure 
distributions  equal  to  the  interfacial  contact  pressures  and  applied 
at  tt  radians  relative  to  the  contact  region,  (see  fig.  6.1).  This 
assumption  is  not  inherent  in  the  method  of  analysis  but  it  does  effect 
the  form  of  the  influence  functions  used  in  the  analysis  to  follow. 
Therefore  the  contact  region  will  be  limited  to  contact  over  half  the 
cylinder,  i.e.,  a half  angle  contact  of  ^ radians  (which  should  cover 
all  cases  of  practical  interest).  A cross  section  of  the  cylinder, 
body  1,  and  cylindrical  seat,  body  2,  is  shown  in  figure  6.2-a. 

The  contour  curves  are  the  same  as 

those  obtained  in  the  analysis  of  the  sphere  and  seat,  therefore  much 
of  the  development  takes  the  same  form.  The  reader  is  referred  to 
chapter  5 for  the  details  of  the  formulation  omi tted  in  this  section. 

It  must  be  remembered  that  the  problem  of  a cylinder  and  seat  in 
contact  is  two  dimensional  and  the  force  applied  at  any  point  on  either 
body  represents  a line  load  with  units  [lbs/inl. 

Surface  points  on  the  cylinder  and  seat  will  be  described  as 
in  the  last  chapter,  i.e.,  coordinates^  pertain  to  points  on  the  cyl- 
inder and  coordinates  0 defines  points  on  the  seat.  The  initial 
location  of  points  A and  B on  the  contour  curves  are  chosen  in  the  same 
way  as  those  on  the  sphere  and  spherical  seat,  s-j  = S2*  The  local 

AAA 

coordinate  system  (r,  t,  <d)  is  fixed  at  point  A as  before,  however,  the 
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Fig. 6. 1.  Cross  section  of  cylinder  and  cylindrical  seat  in 
conformal  contact. 
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Fig.  6.2-a.  (repeated)  Conformal  contact  between  a sphere  and 
spherical  seat  or  a cylinder  and  cylindrical  seat. 
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unit  vector  co  is  parallel  to  the  axis  of  the  cylinder. 

The  contact  criterion  of  equation  (5.3)  exactly  represents 
the  criterion  needed  for  the  solution  of  this  problem,  i.e., 


Sf  - f /C&l  f - &/C44CC  + (W,+Wx)*44  A (5.3) 

where  f is  the  profile  function,  6 is  the  approach  and  w-j  and  w^ 
represent  the  displacements  in  the  radial  directions  of  points  A and 
B on  the  cylinder  and  seat  respectively.  The  quantities  a.  A,  6 
and  £ are  defined  by  equations  (5.1),  (5.4),  (5.5)  and  (5.6)  respec- 
tively. 

A A A A 

Now  consider  the  (n-| , ^ ) and  (r^*  r,^)  cartesian  coordinate 
systems  fixed  to  the  cylinder  and  seat  respectively  as  illustrated  in 
figure  6.2.  Point  A is  located  on  the  cylinder  by  polar  coordinates 
(i|>,  R-j).  The  elastic  displacement  w-j  at  point  A,  can  be  expressed  as 
a function  of  the  pressure  distribution  p (^)  by 


KO,  i'oO  = /?,}  ftriG.tv ; r / m , f, )oiy ' (6-n 

Jl, 

where  G-j  (i|>,  , E-j ) is  given  by  G in  equation  (4.40)  and  ^ is  the 

contact  surface  on  the  cylinder.  Similarly,  on  the  cylindrical  seat, 
the  elastic  displacement  w^  may  be  defined  by 


w2(t>  ) = Rt  f W)  vtA)d*' 


(6.2) 


CYL  INDER 
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CYLINDRICAL  SEAT- 


Fig.  6.2.  Coordinate  systems  for  cylinder  and  cylindrical  seat 
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where  (0,  01  , v^,  E^)  is  defined  by  equation  (4.53)  and  is  the 
contact  region  on  the  cylindrical  seat. 

Combining  equations  (6.1),  (6.2)  and  (5.3)  the  contact 
criterion  may  be  written 


The  suppl ementary  conditions  associated  with  equation  (6.3)  are 
equations  (5.10  a-d) . 


The  initial  choice  of  the  tentative  contact  region  is 
based  on  the  same  assumptions  used  in  the  previous  chapter,  i.e., 
s-|  = s^  so  that 


A simply  discretized  solution  of  equation  (6.3)  is  found  by  first 

★ 

subdividing  fi-j  and  ^ into  a large  number  (N)  of  infinitely  long 
cells  which  are  oriented  such  that  the  infinite  dimensions  of  the 
cells  are  parallel  to  the  axis  of  the  cylinder.  The  normal  pressure 
distribution  is  then  approximated  with  a piecewise  constant  pressure 
distribution  such  that  the  pressure  in  each  cell  is  constant.  The 


(6.3) 


* 


r* MX  “ ft  7/M4X 


(6.4) 
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f.  L. 

iL  cell  on  the  cylinder  will  be  denoted  ft, ..  while  that  on  the  seat 

* 

will  be  termed  ftg..  The  choice  of  is  such  that  it  merges  with 

ftl  ■ . This  is  achieved  by  first  choosing  the  cells  on  ft-j  and  then 

using  the  point  mating  procedure  to  determine  the  cell  boundaries  on 
★ 

ft^  . Cell  i will  be  located  on  the  cylinder  between  ^ and  ip.  + -j 

while  cell  i on  the  seat  will  be  defined  between  <|>.  and  <|>.  + .|.  N 

field  points  are  chosen  on  each  body  such  that  one  lies  within  each 

cell.  The  N + 1 field  point  was  located  on  the  boundary  of  the 

h 

contact  region  on  each  body.  The  location  of  the  l field  point  on 

the  seat  is  determined  using  the  point  mating  procedure  so  that  it  is 

t h 

assumed  to  merge  with  the  i field  point  on  the  cylinder.  Equation 
(6.3)  may  be  written  in  discretized  form  as 


St  = T/COt  j*  ~ S /tOt  Ot  + 


(6.5) 


<5/ 

+ 404  A Z Pt- <)-->, ««) 


The  integrals  in  equation  (6.5)  may  be  evaluated  analytically  as 
described  in  the  following  section.  Thus  N + 1 equations  are  gener- 
ated in  N + 1 unknowns  which  may  be  solved  for  by  Gaussian  elimination 
as  discussed  in  section  5.2. 
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It  remains  to  check  the  validity  of  the  assumption  used  to 
determine  the  location  of  merging  points.  This  is  done  in  precisely 
the  same  manner  as  described  in  section  5.2  using  the  point  mating 
procedure.  The  separation  in  the  tangential  direction  is  determined 
using  equation  (5.13a),  however,  the  tangential  displacements  u-j  and 
U2  used  in  equation  (5.13a)  are  now  defined  by 

u, (vj)  = £ pt  t,J H, ('thy\  (6.6) 

jiu 

and 

Ui  (<Pj)  -•  2 Pi  Plf  Hz(  <f>Jj  <j>\  , Fj  Cty'  (6.7) 

with  and  defined  by  H in  equations  (4.41)  and  (4.54)  respec- 
tively. The  remaining  steps  in  the  point  mating  procedure  are  performed 
as  described  in  section  5.2. 

For  the  two  dimensional  problem  of  plane  strain,  the  only 
meaningful  strains  which  can  be  calculated  (from  the  surface 
displacements)  are  for  the  cylinder  and  e for  the  seat.  The 

strains  e , e , e,  and  e.  are  all  identically  zero  while  e .,  e . 

oxjo  rto  lpo)  4>(jo  n|)  r4> 

and  e can  not  be  calculated  with  only  knowledge  of  the  surface  point 
displacements,  and  e can  be  calculated  using  equations  (5.16) 
and  (5.18)  respectively. 
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6.3  Numerical  Procedures 

In  each  simply  discretized  solution  the  boundaries  on  and 

were  defined  by  ^max  and  <j>  respectively.  The  contact  region  on 

each  body  was  partitioned  into  N cells  by  dividing  a or  ih  by  N 

3 ymax  wnax 

f*  h 

and  defining  the  iL  cell  to  lie  between  ^ and  + -j  where 
^ + 1 - ^ and  *1  + , = i*4max  / N)  (see  fig.  6.3).  The 

pressure  distribution  between  ip  = 0 and  \p  = - ^ was  assumed  to  be 

h 

symmetric  with  respect  to  the  n-|  axis.  The  iLn  field  point  on  the  cyl- 

^i  + ^i  + 1 

inder  was  located  on  the  contour  curve  at  an  angle  of ^ • 

The  N + 1 field  point  was  located  at  + y Similar  locations  for 
the  field  points  on  the  seat  were  chosen  in  terms  of  <J>. , + -j  and 

+ T 

The  integrals  in  equations  (6.5),  (6.6)  and  (6.7)  were 
evaluated  using  the  analytic  formulations  derived  in  appendix  0.  Since 
the  pressure  distribution  is  symmetric  about  the  n-|  axis  (fig.  6.3), 

± . j_ 

the  integral  over  the  i cell  on  the  cylinder  (or  seat)  consisted  of 
two  parts,  the  region  between  i|>.  and  + -j  (or  <|>.  and  <J>.  + i ) and 
between  - ij>.  and  - if;.  + i (or  - <|>.  and  - <!>..  + ^ ) . In  evaluating  the 
tangential  displacements  in  equations  (6.6)  and  (6.7)  the  integrals 
were  multiplied  by  either  + 1 or  - 1 depending  on  the  relative  posi- 
tions of  the  field  point  and  regions  of  integration.  If  the  region  of 
integration  was  located  to  the  left  (in  the  - ip  or  - <p  direction)  of 
the  field  point  the  integral  was  multiplied  by  - 1 ; otherwise  it  was 
multiplied  by  1.  This  procedure  accounted  for  the  sign  of  the 
direction  of  the  displacement  due  to  the  position  of  the  loading.  The 
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integrals  when  evaluated  at  their  singularities  present  no  problem 
since  the  results  arefinite  as  shown  in  appendix  0. 

The  point  mating  procedure  was  employed  in  the  same  manner  as 
described  in  section  5.3,  to  converge  on  the  coordinates  of  a set  of 
merging  field  points.  Also  as  described  in  section  5.3,  N + 1 
equations  were  generated  and  the  variables  were  solved  for  using 
Gaussian  elimination. 

The  total  applied  force  per  unit  length  can  be  calculated 
from  the  discretized  pressure  distribution  determined  in  the  analysis. 
Consider  the  constant  normal  pressure  P.  over  cell  i on  the  cylinder. 
The  component  of  incremental  force  in  the  - n-j direction  at  angle  ip 
is  P.j  cos  ip  R-j  dip  . Integrating  this  between  ^ and  ^ + -j  the  force 

per  length  over  cell  i becomes  P^  (sin  ip . + -|  - sinij;  .).  Recalling 
that  the  pressure  acts  over  an  identical  region  between  - ijj.  and 
- ijj.  + i»  the  total  force  on  the  cylinder  may  be  found  by  summing  the 
forces  on  each  cell,  i.e., 

N 

P-  2RS  Jf  p-  ( s&on  Yc  ) (6-8) 

C'  i 

Finally  the  approach  6 may  be  calculated  by  back  substitution  of  the 
P , s in  the  N + Inequation  written  at  the  N + 1 field  point  (see 
chapter  2 for  a description  of  the  complete  simply  discretized  method 
of  solution). 
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6.4  Numerical  Results 

The  following  numerical  example  in  plane  stress  was 
considered: 

R.j  = 1.00  in 

R2  = 1.01  in 

"-max  = 40  Deg  (6.9) 

E]  * E2  * 30  x 106  psi 

Vi  = v2  = 0.3 

A program  C0NCYL  was  written  using  the  analysis  of  section 
6.3.  The  results  for  the  above  problem  were  compared  to  those  of 
Persson  [19641  . 

The  pressure  distributions  obtained  by  Persson  and  C0NCYL 
are  plotted  in  figure  6.4.  A close  correspondence  exists  between  the 
two  solutions.  The  displacements  calculated  in  C0NCYL  are  tabulated 
in  tables  6.1  and  6.2  as  functions  of  the  angles  ip  and  <J>. 
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Fig.  6.4.  Pressure  distribution  between  a cylinder  and  cylindrical  seat,  ^max 
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TABLE  6.1 

DISPLACEMENTS  ON  CYLINDER 


[Deg] 

Radial  3 

Displacement  * 10J 
[in] 

Tangential  ^ 

Displacement  * 10J 
[in] 

1.33 

2.019 

0.1101 

4.00 

2.006 

0.3288 

6.67 

1.979 

0.5429 

9.33 

1.939 

0.7492 

12.00 

1.886 

0.9450 

14.67 

1.819 

1.127 

17.33 

1.739 

1.294 

20.00 

1.646 

1.443 

22.67 

1.540 

1.573 

25.33 

1.422 

1.682 

28.00 

1.290 

1.771 

30.67 

1.146 

1.841 

33.33 

0.9885 

1.896 

36.00 

0.8189 

1.945 

38.67 

0.6368 

2.010 

40.00 

0.5409 

2.039 

PLANE  STRESS  RESULTS  FOR: 

^max  = 40°  * R1  = 1*00  in’  R2  = 101  in’ 
E-|  = E^  = 30  x 10^  psi,  v.j  = = 0.3 
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TABLE  6.2 

DISPLACEMENTS  ON  CYLINDRICAL  SEAT 


<p 

[Deg] 

Radial  3 

Displacement  * 10J 
[in] 

Tangential  - 

Displacement  * 10J 
[in] 

1.31 

0.3097 

-0.07100 

3.93 

0.3078 

-0.2124 

6.55 

0.3040 

-0.3521 

9.17 

0.2984 

-0.4889 

11.79 

0.2910 

-0.6216 

14.41 

0.2817 

-0.7490 

17.04 

0.2706 

-0.8700 

19.66 

0.2577 

-0.9834 

22.28 

0.2431 

-1.088 

24.91 

0.2268 

-1.182 

27.54 

0.2089 

-1.265 

30.16 

0.1894 

-1.334 

32.79 

0.1683 

-1.389 

35.43 

0.1458 

-1.427 

38.06 

0.1219 

-1.444 

39.38 

0.1094 

-1.446 

PLANE  STRESS  RESULTS  FOR: 

= 40°  , R,  = 1.00  in,  R„  = 1.01  in, 

Mid  X S c. 

E1  = E9  = 30  x 106  psi,  v = v2  = 0.3 
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The  radial  displacements  on  the  seat  are  larger  by  a factor 
of  1.5  to  2.0  than  those  on  the  cylinder.  As  in  the  case  of  the  sphere 
and  seat,  the  tangential  displacements  on  the  cylindrical  seat  were 
found  to  be  negative  while  the  tangential  displacements  on  the  cylinder 
were  positive.  The  physical  interpretation  of  this  result  is  the  same 
as  that  expressed  in  section  5.4  for  the  sphere  and  seat.  The  values 
of  force  and  approach  for  this  problem  were  found  to  be 
F = 0.6922  x io5  lbs/in  and  6 = 0.5120  x 10~2  in. 

In  order  to  correlate  results  with  Hertzian  theory,  the 
problem  for  = 0.1  degree  was  analyzed.  The  resulting  pressure 
distribution  is  plotted  in  figure  6.5  along  with  the  results  of  Persson 
and  Hertz.  There  is  close  agreement  between  all  solutions  as  would  be 
expected  for  this  case  of  small  contact  area. 

Figure  6.6  illustrates  the  relation  between  the  load  F,  radial 

difference  ArJ  and  the  maximum  angle  of  contact,  ip  . Along  with 

md  x 

the  results  of  C0NCYL  are  plotted  the  solutions  of  Hertz,  Sjtaerman 
[1949]  and  Persson. 

There  is  close  agreement  between  all  solutions  for  angles  of 
contact  less  than  15  degrees.  For  larger  angles  of  contact  there  is  a 
close  correspondence  between  the  results  of  C0NCYL  and  those  of  Persson. 
The  curves  corresponding  to  Hertz's  theory  and  that  of  Sjtaerman 
deviate  significantly. 

Figure  6.7  illustrates  the  variation  of  maximum  pressure  with 

^ AR  equals  R^  - R^ . 


(-t-)J  ‘X 
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Fig.  6.5.  Pressure  distribution  between  a cylinder  and  cylindrical  seat, 
0.100. 
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R p ^ 

Fig.  6.7.  Relationship  between  HM  and  tangent  of 
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maximum  contact  angle  for  both  the  present  theory  and  that  of  Persson. 
Figure  6.8  relates  the  variation  between  the  maximum  pressure,  radial 
difference  and  load.  In  each  of  the  last  two  figures  the  results  of 
CONCYL  show  close  agreement  with  the  results  of  Persson. 

The  computer  costs  in  running  CONCYL  were  minimal,  being  about 
twenty-three  cents  for  a fifteen  node  case.  This  corresponds  to  about 
five  seconds  of  CPU  time  on  the  IBM  360/65  computer.  The  low  costs  in 
CONCYL  can  be  attributed  to  the  fact  that  all  integration  was  performed 
analytically  rather  than  numerically  as  in  CONSPHERE. 

6.5  Conclusions 

It  can  be  concluded  that  the  problem  of  a cylindrical  seat 
has  been  successfully  solved  using  the  conformal  theory  presented  in 
chapter  3.  The  pressure  distributions  for  the  problems  where 
w = 40°  and  = 0.1°  were  found  to  agree  with  the  results  of 
Persson  [1964],  In  addition  the  latter  results  also  corresponded  to 
the  solution  of  Hertzian  theory.  The  displacements  were  calculated 
for  the  case  where  \pmQx  = 40°.  The  tangential  components  were  found 
to  be  of  opposite  sign  on  each  body.  The  load  vs.  subtended  angle 
relationship  was  found  to  agree  with  the  solution  of  Persson,  however, 
the  solution  of  Sjtaerman  [1949]  deviates  significantly  for  angles 
greater  than  20  degrees.  The  results  of  maximum  pressure  vs.  contact 
angle  also  agreed  well  with  that  of  Persson. 

In  general  there  was  strong  agreement  with  the  Persson  solu- 
tion. This  close  agreement  supports  Persson's  assumption  that  the 
contour  curve  of  the  contact  region  is  circular.  The  correlation 
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R P 

Fig.  6.8.  Relationship  between  ^ max  and 
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between  the  solutions  presented  and  the  Hertzian  solutions  broke  down 
for  large  angles  as  expected.  The  results  of  Sjtaerman  also  deviated 
from  those  of  the  present  theory  which  may  be  explained  by  the  fact 
that  Sjtaerman  may  have  used  an  extremely  crude  finite  difference 
scheme  since  he  could  not  take  advantage  of  modern  computational  aids 
at  the  time  of  publication. 


7.  CONTACT  STRESSES  FOR  MULTIPLY  CONNECTED  REGIONS 


7.1  Introduction 

Contact  problems  involving  multiply-connected  contact 
regions  have  received  little  attention  in  the  literature,  possibly 
because  of  the  non-Hertzian  nature  of  such  problems.  Such  problems 
arise,  for  example,  whenever  either  of  the  contacting  bodies  have 
surface  pits  (e.g.,  casting  defects,  corrosion  pits,  machining 
faults,  etc.).  Barely  perceptible  surface  flaws  can  cause  high 
stress  concentrations,  and  consequently,  rapid  fatigue  failure. 
Experimental  observations  by  Tallian  [1967],  Martin  and  Eberhardt 
[1967]  and  Littman  and  Widner  [1966]  indicate  that  such  surface 
defects  may  be  potential  nuclei  of  microcrack  propagation  and  can 
produce  rapid  destruction  of  rolling  surfaces. 

Based  on  the  degree  of  difficulty  associated  with  their 
solution,  these  problems  may  be  divided  into  the  following  two 
categories: 

(i)  Contact  region  known  a priori: 

When  the  indentor  contact  surface  is  flat  (or  almost  flat)  it 
will  be  called  a "stamp,"  and  the  contact  surface  is  defined 
a priori  by  the  stamp  boundary.  When  the  indentor  surface 
is  not  flat,  but  the  indentor  has  a substantially  higher 
elastic  modulus  than  the  indented  body,  the  indentor  can  be 
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treated  as  rigid,  and  the  shape  of  the  contact  region  becomes 
known  for  any  given  depth  of  penetration  relative  to  the 

indentor  tip.  These  are  also  termed  "punch"  problems. 

(ii)  Elastic  Contact  Problems: 

When  the  indentor  is  not  a stamp,  and  the  two  bodies  have 
comparable  elastic  moduli,  then  the  geometry  of  the  contact 
region  is  unknown  a priori,  and  it  must  be  determined  by 
solving  the  appropriate  Elasticity  problem. 

To  the  best  of  our  knowledge,  only  one  recent  solution  by 
Chaud  et  al.  [1974]  for  three  dimensional  elastostatics  with  multiply- 
connected  regions,  has  been  reported  in  the  literature.  However, 
solution  of  a few  special  cases  of  rigid  indentor  problems  (category 
[i])  have  been  found  by  Olesiak  [1965],  Parlas  and  Michalpoulos  [1972] 
and  Chiu  [1969]. 

Olesiak  [1965]  solved  the  problem  of  an  annular  flat  faced- 
stamp  pressed  on  an  elastic  half  space.  Parlas  et  al . proposed  the 
solution  for  a "bolt  shaped"  indentor  pressed  into  an  elastic  half 
space  with  a cylindrical  hole.  The  cylindrical  (bolt)  section  of 
the  indentor  was  assumed  to  be  rigidly  bonded  to  the  wall  of  the 
cylindrical  hole  while  the  bottom  face  of  the  bolt  head  presses 
against  the  half  space. 

Chiu  [1969]  solved  the  problem  of  an  infinitely  long  rigid 
cylinder  in  contact  with  an  elastic  half  space,  where  the  rigid 
cylinder  has  a groove  running  parallel  to  its  axis. 

In  this  chapter,  results  indicate  that  problems  of 
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both  categories  (i)  and  (ii)  may  be  successfully  solved  by  an  extension 
of  the  method  introduced  by  Singh  and  Paul  [1974]. 

A brief  synopsis  of  the  formulation  and  application  of  the 
"simply  discretized"  method  of  solution  are  given  with  some  limitations 
and  advantages  of  this  method  in  section  7.2.  The  example  problem  of  a 
pitted  sphere  in  contact  with  a complete  sphere  is  described  in  section 
7.3.  Techniques  devised  for  an  accurate  numerical  solution  and  rapid 
convergence  are  described  in  section  7.4.  Results  for  an  example  are 
given  in  section  7.5,  and  conclusions  are  reviewed  in  section  7.6. 

7.2  Formulation 

We  will  restrict  our  attention  to  "nonconformal " contact 
problems  where  the  dimensions  of  the  contact  region  are  small  compared 
to  appropriate  radii  of  curvature  of  the  undeformed  bodies.  There- 
fore,  we  may  assume  that  the  contact  surfaces  do  not  deviate  signifi- 
cantly from  a reference  plane  in  which  we  imbed  fixed  cartesian 
axes  (x,y).  Furthermore,  we  shall  consider  only  those  cases  where 
the  two  bodies  undergo  a relative  rigid  body  translation  of  amount 
6,  in  a direction  normal  to  the  reference  plane,  plus  an  elastic 
deformation.  The  translation  6 is  called  the  "relative  approach"  and 
is  positive  if  it  moves  the  bodies  towards  one  another.  We  will  also 
assume  that  the  applied  load  consists  of  a force  F,  acting  normal 
to  the  reference  plane,  and  that  the  contacting  surfaces  have  a 
sufficient  degree  of  symmetry  that  the  resultant  of  the  contact  pres- 
sures on  each  body  is  a force  of  magnitude  F which  acts  through  the 
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origin  0 of  the  reference  plane  and  equilibrates  the  applied  force  F. 

The  fundamental  integral  equation  governing  nonconformal 
contact  problems  was  shown  in  chapter  2 to  be 


St*.?) 


C P(K*')JtW  x r 

JSL 


(7.1) 


where  the  "elastic  parameter"  k is  defined  as 


(7.2) 


In  the  foregoing  equations,  and  denote  the  Poisson's 

ratio  and  Young's  modulus  respectively  for  body  1 (indentor)  and  body 
2 (indented);  p(x',y')  is  the  normal  pressure  over  the  contact 
surface;  ft  is  the  projection  of  the  contact  surface  on  the  (x,y) 
reference  plane;  f(x,y)  represents  the  initial  separation  (or  gap) 
between  surface  points  on  the  two  bodies,  located  at  the  same  (x,y) 
coordinates,  before  the  load  F is  applied;  S(x,y)  is  the  separation 
of  the  opposed  surface  points  after  the  load  is  applied.  Figure  7.1 
illustrates  the  initial  separation  f for  a case  of  axial  symmetry 
where  f is  a function  f(r)  of  the  radial  coordinate  r. 

The  condition  of  impenetrability  of  matter  requires  that 
S(x,y)  should  vanish  inside  ft  and  it  should  be  positive  outside  of 
ft.  Conversely,  the  interfacial  contact  pressure  p(x,y)  should  be 
positive  inside  ft,  and  it  should  vanish  identically  outside  of  it. 
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In  symbolic  terms. 


S'O  foR  (X,  V)  tNStpe  AhPOAiA.  (7.3a) 


S ?0  toR(x,y)  ourstpe  ji 


(7.3b) 


P(^t^)  sO  fOR(X,V)  OVTStDEA  (7.4a) 

PCKjV)  ZO  totCx,*)  iN$/pe /y*poMSL  (7.4b) 


(7.4a) 


In  short,  a solution  of  the  problem  requires  the  determina- 


tion of  the  boundaries  of  region  ft,  a pressure  field  p(x,y),  and  an 
approach  6 which  satisfy  relations  (7.1 )-(7.4).  The  associated 
load  may  be  found  from  the  expression 


The  absence  of  foreknowledge  of  the  contact  region  ^ is  a major 
impediment  to  a mathematical  solution.  This  obstacle  is  overcome  by 
postulating  a tentative  contact  region  ft*.  Singh  and  Paul  [1974] 
proposed  that  the  "interpenetration  curve"  described  by 


be  used  as  a tentative  contact  region.  Equation  (7.6)  defined  the 
contour  of  the  curve  formed  by  interpenetration  (without  deformation) 
of  the  two  surfaces  through  an  arbitrary  distance  d.  Picking  a 
suitable  value  of  d establishes  the  candidate  contact  region  ft*. 
Using  this  as  a preliminary  estimate  of  ft,  equation  (7.1)  is  readily 


A 


(7.5) 


(7.6) 


recognized  to  be  an  integral  equation  of  the  first  kind. 


Equation  (7.1)  can  be  solved  using  the  "simply  discretized" 
method  of  Singh  and  Paul  which  is  reviewed  in  chapter  2.  A "simply 
discretized"  numerical  solution  of  equation  (7.1)  is  found  by 
subdividing  ft  into  a large  number  of  small  cells.  The  pressure 
function  p(x,y)  is  replaced  by  a piecewise  constant  pressure  field 
(pressure  P.  in  cell  i).  Thus  if  ft  is  subdivided  into  N cells, 
equation  (7.1)  becomes 


where  ft..  is  the  region  of  cell  i.  In  equation  (7.7),  N values  of  p^ 
and  the  constant  6 are  unknowns  to  be  determined.  The  centroids 
(x. ,y. ) of  the  cells  are  taken  as  field  points  (x,y)  and  equation 
(7.7)  is  written  for  each  field  point.  The  integrals  in  equation 
(7.7)  are  evaluated  by  numerical  quadrature.  Thus  N linear  algebrai 
equations  are  generated.  An  additional  independent  linear  equation, 
essential  for  a unique  solution,  is  generated  by  picking  up  a field 
point  other  than  the  cell  centroids.  The  choice  of  this  additional 
field  point  is  otherwise  arbitrary,  however,  it  does  affect  the 
quality  of  the  results,  as  discussed  in  section  7.4. 

Having  thus  generated  a set  of  N + 1 linear  equations,  the 
N unknown  pressures,  P.,  and  the  approach  6,  are  obtained  through 


(7.7) 
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Gaussian  elimination.  The  next  step  in  the  solution  is  to  determine 

★ 

whether  the  tentatively  selected  region  of  integration  ft  is  indeed 
the  true  contact  region.  This  is  done  by  utilizing  the  inequalities 
(7.3)  and  (7.4)  and  systematically  adjusting  the  boundaries  of  ft 
until  these  inequalities  are  satisfied. 

Singh  and  Paul  [1974]  showed  that  the  "simply  discretized" 
method  was  unstable  in  the  general  case  and  was  incapable  of  pre- 
dicting the  proper  stress  distribution.  For  such  problems  they  found 
it  necessary  to  introduce  stabilizing  techniques  known  as  the 
"Redundant  Field  Point  Method,"  and  the  "Functional  Regularization 
Method"  (see  Singh  [1972],  Singh  and  Paul  [1973]). 

The  amount  of  numerical  computation  required  for  either  of 
the  two  last  methods  exceeds  that  of  the  Simply  Discretized  Method. 
Accordingly,  it  is  desirable  to  use  the  latter  whenever  circumstances 
permit. 

In  this  chapter  we  will  focus  on  a problem  with  complete 
axisymmetry,  and  it  will  be  shown  that  the  Simply  Discretized  Method 
provides  an  excellent  solution,  provided  that  the  maximum  possible 
use  is  made  of  the  symmetry  of  the  problem. 

In  other  words,  we  recognize  that  all  cells  located  at  the 
same  radius  from  the  axis  of  symmetry  have  the  same  contact  pressure 
at  their  centroids,  and  the  number  of  unknown  pressures  P^  is  reduced 
from  the  number  of  cells  to  N (the  number  of  annular  rings  formed 
by  an  axisymmetric  distribution  of  cells).  By  using  the  Simply 
Discretized  Method,  we  are  able  to  utilize  Inequality  (7.4)  to 
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iteratively  refine  the  region  of  contact  Upon  satisfying 
Inequality  (7.4),  it  was  invariably  found  that  Inequality  (7.3)  was 
satisfied. 

The  nature  of  the  Functional  Regularization  Method  prohibits 
the  use  of  Inequality  (7.4)  as  a basis  for  refining  n. 

Numerical  experiments  have  indicated  that  iteration  pro- 
cedures based  on  Inequality  (7.4)  converge  much  faster  than  those 
based  upon  Inequality  (7.3).  Further  details  of  the  iteration  pro- 
cedures will  be  found  in  sections  7.4  and  7.5. 

7.3  Pitted  Sphere  Geometry 

As  a typical  example,  contact  of  a pitted  elastic  sphere  of 

radius  with  an  unpitted  elastic  sphere  of  radius  R^  is  considered. 

A section  of  the  pitted  surface  by  a plane  through  the  axis  of 

symmetry  is  shown  in  figure  7.1.  The  local  contour  of  the  pitted 

surface  is  idealized  as  a torus  smoothly  blended  into  a sphere. 

The  blending  point  P,  where  the  pit  joins  the  main  surface,  is 

located  at  a distance  r,  from  the  load  line.  The  center  of  curvature 

b 

I 

0 of  the  pit  blending  arc  lies  on  the  conical  surface  of  semi  vertex 

angle  x.  The  meridional  radius  of  curvature  of  the  torus  is  r . 

Note  that  the  discontinuity  in  curvature  which  occurs  at  P 
does  not  preclude  the  use  of  the  method  of  solution  being  used. 

A tentative  contact  region,  ft,  is  established  by  a hypothetical  inter- 
penetration of  the  two  spheres  through  a distance  d.  The  annulus  of 
contact  so  formed  is  bounded  by  an  inner  radius  r^  and  an  outer  radius 
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Fig.  7.1.  Geometry  of  pitted  surface 


i 


Fig.  7.2.  Generation  of  annular  interpenetration 
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r as  shown  in  figure  7.2,  where  suitable  coordinate  axes  r,  and  z 
o 

are  indicated.  The  values  of  r^  and  rQ  for  a given  problem  are 
determined  as  follows.  The  z coordinate  of  a point  C(p,z-j)  located 
at  a distance  p from  the  z-axis  on  the  toroidal  portion  of  body  1 
(see  fig.  712),  where 


/<  rk 


(7.8a) 


is : 


z.  = t, -fa-*)  r-V'TwT-T)1 


where 


(R,-re)H 

R, 


■ -i 


(7.8b) 


(7.8c) 


(7.9) 


The  z-coordinate  of  a point  on  sphere  2,  located  at  a distance  p from 
the  z-axis  is  given  by 


(7  io) 

Since  point  C lies  on  both  the  torus  and  the  lower  sphere,  z-j  = 
thus  equations  (7.8b)  and  (7.10)  require  that 
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cl  - R)  - (R,~rc)  CfrL  T 

P<rw  <7-lla> 

Furthermore,  the  z-coordinate  of  a material  point  C'  located  on  the 
spherical  portion  of  body  1 , at  a distance  p from  the  z-axis,  is  given 
by 


* *>'  ’S')  1 (7.12) 

where 

P > 1%  (7.12a) 

Hence,  for  a given  interpenetration  d,  the  radius  p of  a point  on  the 
intersection  of  sphere  2 and  spherical  region  of  body  1 is  given  by 

d *•  R,  - (K, -(e'-f)  'H  (7.1 3a ) 

P>t%  (7.13b) 

The  geometry  of  the  toroidal  surface  indicates  that  for  rc  < R] > 
equation  (7.11)  has  two  solutions  for  p.  Let  p-j  and  p2  (p-|  < p^)  be 
roots  of  equation  (7.11).  Two  cases  are  readily  identified. 

Case  (i).  When  both  inner  and  outer  radii  of  the  assumed 
contact  region  lie  inside  the  blending  radius,  i.e., 

A <n 


(7.14) 
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In  this  case  the  contact  is  assumed  to  be  completely  confined  to 
the  toroidal  segment  of  body  1,  in  which  case 

~ P\ 

(7.15) 

*0  r ft 

Case  (ii).  When  the  outer  boundary  of  ft  lies  beyond  the 
blending  radius  (as  shown  in  fig.  7.2),  i.e., 

P7>h  (7.16a) 

In  this  case 


(7.16b) 


and  the  outer  radius  rQ  is  determined  from  solution  of  equation  (7.13). 
Note  that  equations  (7.11)  and  (7.13)  are  transcendental  in  p,  which 
can  be  found  by  an  iterative  procedure  (e.g.,  Newton  Raphson). 

In  order  to  find  the  initial  separation  f(r),  shown  in  figure 
7.1,  it  is  only  necessary  to  find 

{(r)  - 2,  - 2Z  (7.i6c) 


where  Z2  is  found  from  equation  (7.10)  with  p = r and  d = dQ;  dQ  is 
the  value  of  d corresponding  to  initial  contact  as  shown  in  figure  7.1. 
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To  find  z.j , set  p,  e r and  use  equation  (7.8b)  for  points  on  the  torus 
(r  < rb),  or  equation  (7.12)  for  points  on  the  upper  sphere  (r  > r^). 

In  order  to  find  the  initial  separation  dQ,  it  is  necessary 
to  note  from  figure  7.1,  that  when  d = dQ,  the  slope  of  the  torus 
matches  that  of  the  lower  sphere  at  the  contact  point;  i.e., 

oil,  _ el 

O IP  dfi  (7.17) 

where  the  derivatives  are  found  from  equation  (7.8b)  and  equation 
(7.10).  Equation  (7.17),  together  with  equations  (7.8b)  and  (7.10), 
suffice  to  find  dQ,  and  the  two  coordinates  (r,z)  of  the  initial 
contact  point. 

Having  found  the  boundaries  (r  and  rQ)  of  the  contact  region 
ft  and  the  initial  separation  function  f(r),  we  may  proceed  to  solve 
the  governing  integral  equation  (7.1). 

7.4  Numerical  Solution  Procedure 

The  contact  region  ft  is  subdivided  into  N annular  rings. 

Since  a steep  pressure  gradient  is  expected  near  the  pit,  the  annular 
rings  near  the  inner  boundary  are  very  narrow  in  width.  It  was  also 
learned  from  experience  that  the  peak  pressure  always  occurs  at 
some  radius  r where  r < r^.  Guided  by  this  consideration,  a majority 
of  the  rings  are  clustered  in  the  region  rj  £ r <_  r^.  Exploiting  the 
axisymmetry  of  the  problem,  we  assume  that  the  pressure  is  constant 
in  each  ring.  The  rings  are  numbered  sequentially  from  1 to  N,  from 
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the  inside  out,  and  the  pressure  in  the  i-th  ring  is  assumed  to  be 

an  unknown  constant  . Let  r.  and  r^+^  be  the  inner  and  outer 

bounding  radii  for  cell  i;  thus  r^  e r^  and  r^  e r . Each  ring 

is  further  subdivided  circumferentially  into  m equal  sectors  by 

drawing  (m)  equispaced  radial  rays  from  the  center  of  fi;  the  angle  A(p 

between  two  adjacent  rays  is  2n/m.  The  sector,  bounded  by  radial 

rays  1 and  2,  is  shown  in  figure  7.3. 

The  region  of  the  sector  located  in  the  i-th  ring,  between 

ray  j and  ray  (j+1),  is  identified  as  S..;  and  its  centroidal  radius 

J 

by  3 — . Elementary  calculations  show  that 


3 * <t> 


(7.18) 


The  centroids  of  the  first  sector  shown  in  figure  7.3  (i.e.,  where 
j = 1)  are  selected  as  field  points. 

Thus  for  the  field  point  H,  equations  (7.7)  and  (7.3a)  reduce 
to 


N it) 

It  £ F 

1*1  J*/  J 


c/Aj 


3; 


iji 


u 


-S 


(7.19) 
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where  f(r)  is  calculated  for  r = 3^  from  equation  (7.16c).  c^ 

is  the  radial  distance  from  field  point  £ to  the  elemental  area  dA.. 

■ J 

located  in  S...  For  most  cells,  the  integral  in  equation  (7.19)  may 
^ J 

be  replaced  by  the  approximation 


«,  - j 


cM 


IJ 


S .. 


» • j 

tji 


Ai 


(7.20) 


where  c.  . is  the  distance  between  field  point  £ and  the  centroid  of 

* J 

the  region  S. .,  whose  area  is  denoted  by  A. ..  It  was  shown  in 
Singh  [1972]  that,  in  general,  equation  (7.20)  is  a very  useful 
approximation  which  results  in  a significant  reduction  of  computa- 
tion time,  without  compromising  the  accuracy  of  results.  However, 
for  regions  located  in  the  immediate  vicinity  of  the  field  point  £, 
the  errors  due  to  the  approximation  (7.20)  may  be  unacceptable. 

To  avoid  such  errors,  I...  is  evaluated  by  numerical  quadrature  within 

1 J 36 

cells  located  near  the  field  point.  The  criterion  which  must  be 
satisfied  in  order  to  use  equation  (7.20)  is 


CiJt  > motx  (rA<t>,  Ar) 


(7.21) 


In  equation(7.21 ) , rA<p  and  Ar  are  the  side  lengths  of  a typical  cell. 
Notice  that  when  the  field  point  £ lies  inside  the  region  S..  . (i.e., 
j = 1,  i = £),  c. . = 0,  and  hence  the  integrand  in  equation  (7.15) 

I J 36 

has  a singularity.  However,  for  such  cases,  an  approximate  analytical 
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solution  for  the  integral  is  readily  constructed  as  shown  in 
Appendix  G. 

In  this  manner,  N linear  equations  corresponding  to  the  N 
field  points  are  generated.  An  additional  linearly  independent 
equation  is  generated  by  selecting  point  P'  at  the  outermost  boundary 
of  the  contact  region  as  field  point  (N  + 1).  The  location  of  this 
additional  field  point  has  a pronounced  affect  on  the  solution,  which 
deteriorates  as  P*  is  moved  inside  the  boundary.  It  is  plausible 
to  assume  that  this  behavior  is  due  to  the  gradual  increase  in  cell 
width  Ar  with  r (see  fig.  7.3)  which  was  introduced  to  keep  the  aspect 
ratio  of  the  cells  from  becoming  excessive.  With  the  cells  so 
designed,  the  location  of  P'  shown  in  figure  7.3  maximizes  the  distance 
between  P'  and  its  nearest  neighboring  field  point.  This  in  turn 
tends  to  maximize  the  amount  of  independent  information  supplied  by 
the  equation  written  for  field  point  P',  and  should  tend  to  minimize 
ill-conditioning  effects  on  the  coefficient  matrix  generated. 

Thus  (N  + 1)  equations  in  (N  + 1)  unknowns  are  generated,  and 
equation  (7.19)  assumes  the  form 


(7.22) 


and  the  equation  using  P'  as  a field  point  becomes 


(7.23) 
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where  f..  is  the  value  of  the  "initial  separation"  function  f(r)  at 
the  field  point  i.  f^  is  the  value  of  f(r)  at  P ' ; and  summation 
from  1 to  N is  henceforth  implied  over  repeated  subscripts.  From 
equations  (7.22)  and  (7.23),  6 may  be  eliminated  to  yield 

Bij  Pj  ' f;  (7-24) 


where 


and 


(7.25) 


(7.26) 


When  equation  (7.24)  is  solved  for  P..,  using  Gaussian  elimination,  the 
resulting  pressure  distribution  is  usually  found  to  predict  negative 
contact  pressures  in  the  immediate  vicinity  of  the  inside  boundary, 
r = rj.  The  axisymmetry  of  the  problems  enables  us  to  maintain  the 
outside  boundary  fixed,  and  iterate  on  the  inside  boundary  where  the 
predicted  pressure  is  incorrect.  The  iteration  scheme  is  best 
explained  with  the  aid  of  the  numerical  example  given  in  section  7.5. 

7.5  A Numerical  Example 

The  following  example  problem  was  considered. 

R1  = R2  = 1 in 

V1  = v2  = 0*3 

E]  = E2  = 30  X 106  lb/in2 
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r = .006  in 
c 

r,  = .00025  in 
b 


The  results  are  presented  in  dimensionless 

form.  Let 

„ 2 R,  Rt 

( 7.27) 

Then,  we  define 

★ 

Dimensionless  pressure  in  ring  i,  P.  = 

k P.  (7.28) 

* kF 

Dimensionless  load,  F = 

R2 

m 

(7.29) 

★ 

Dimensionless  distance  from  origin  of  ft,  r = r/R  (7.30a) 

3 m 

'fc 

Dimensionless  approach,  6 = 6/Rm  (7.30b) 


r.  = r,  /R 
b b m 

(7.31a) 

r = r /R 
c c m 

(7.31b) 

Figure  7.4  shows  the  pressure  distribution  near  the  inside  boundary 
for  the  uniterated  solution.  The  pressure  distribution  far  from  the 
pit  agrees  closely  with  the  Hertzian  solution  for  unpitted  spheres 
(not  shown  in  the  figure).  However,  the  pressure  in  cell  #1  is  highly 
negative.  The  pressures  in  the  successive  cells  are  less  and  less 
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negative,  until  at  point  the  pressure  curve  crosses  the  r*-axis. 

The  shape  of  the  pressure  curve  readily  suggests  the  iteration  scheme. 
The  new  region  of  integration  is  assumed  to  have  inner  radius 
rj  = 0Q.| . The  discretized  equation  set  (7.24)  is  generated  corre- 
sponding to  this  new  region  and  thus  a new  pressure  vector  is 
generated  (see  first  iteration,  fig.  7.4).  This  new  curve  also  has  a 
negative  peak  (weaker  than  that  of  the  uniterated  solution)  at  the 
innermost  field  point.  The  new  point  of  intersection  is  Q^,  which 
defines  the  inner  boundary  of  fi  for  the  next  iteration.  The  process 
is  thus  continued  until  all  pressures  are  positive.  In  figure  7.4, 
the  third  iteration  yields  the  desired  solution.  It  is  found  that 
this  solution  also  satisfies  Inequality  (7.3),  thus  qualifying  as  the 
"true"  solution  of  the  contact  problem.  The  complete  pressure  distri- 
bution is  shown  in  figure  7.5.  Notice  the  essentially  Hertzian 
pressure  distribution  (corresponding  to  contact  of  unpitted  spheres) 
at  r*  > 6 x 10’4.  Thus  the  effect  of  the  cavity  is  of  a strictly 
localized  nature.  However,  as  the  cavity  is  made  larger  (e.g., 
rj/rQ  0.3)  the  pressure  curve  departs  completely  from  the  Hertzian 
case.  For  example,  figure  7.6  shows  a typical  pressure  distribution 
for  r j / rQ  = 0.623,  along  with  the  Hertzian  solution  for  unpitted 
spheres  corresponding  to  identical  values  of  thrust  F. 

In  order  to  establish  confidence  in  the  solution,  it  is 
necessary  to  study  its  convergence  with  change  in  the  number  of  cells 
used.  It  must  be  recognized  that  it  is  necessary  for  the  cells  to 
be  densely  concentrated  only  in  that  region  where  a high  pressure 
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Fig.  7.4.  Boundary  iteration  sequence 
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gradient  exists.  Therefore,  for  purposes  of  convergence  studies, 

we  have  systematically  varied  the  number  of  cells  within  a fixed 

radius  r . This  radius  is  chosen  arbitrarily  for  each  problem  in 
con 

such  a way  that  the  major  area  of  stress  concentration  lies  inside 
the  radius  rcor).  For  the  example  problem  considered,  rcQn  = .0003. 

Let  N „ be  the  number  of  rings  located  within  radius  r . Figure 

7.7  illustrates  the  convergence  of  the  peak  pressure,  p*ax-  Figure 

7.8  shows  the  convergence  of  stress  concentration  factor  with  N . 

Stress  Concentration  Factor  (SCF)  is  defined  as  the  ratio  of  the  peak 

computed  pressure  to  the  peak  pressure  for  unpitted  spheres  under 

equal  thrust.  Notice  both  figures  7. 7 and  7. 8exhibi  t convergence  for 

ILn  > 8. 
con 

The  load-approach  curve  is  shown  in  figure  7.9.  It  is  obvious 
from  figure  7.9  that  the  compliance  characteristics  of  the  balls 
(with  small  pits)  remain  essentially  the  same  as  that  predicted  by 
the  Hertzian  solution. 

Figure  7.10  shows  SCF  as  a function  of  cavity  edge  radius 

r*.  Smaller  values  of  r*  cause  greater  stress  concentration.  Due  to 
c c 

the  nonlinearity  of  the  problem,  the  SCF  is  also  a function  of  the 
applied  load  F*.  Table  7.1  shows  the  variation  of  SCF  with  the  size 
of  the  pit  (measured  by  blend  point  radius).  Notice  that  SCF  increases 
with  increasing  value  of  r£.  This  variation  of  SCF  with  r£  may  be 
related  to  the  loss  of  load  carrying  area. 

The  computer  program  developed  to  solve  this  problem  is 
moderately  efficient.  For  example,  the  nine  cases,  needed  to  generate 
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TABLE  7.1 

DEPENDENCE  OF  STRESS  CONCENTRATION  FACTOR  ON 


Case  No. 

r*  x lo3 

D 

SCF 

F*  x IQ8 

o 

X 

r*  x 103 

1 

0.25 

1.692 

0.9743 

0.1023 

0.1845 

2 

0.35 

1.856 

0.9737 

0.1029 

0.2753 

3 

0.50 

2.049 

0.9702 

0.1041 

0.4166 

= r. 

r*  = 0.006,  r*  = 
c 0 

0.002236, 

E = 30  x 106 

psi , v =0.3. 
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Fig.  7.7.  Convergence  of  peak  pressure  with  increasing 
number  of  cells. 
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Fig.  7.8.  Convergence  of  stress  concentration  with 
increasing  number  of  cells. 
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Fig.  7.10.  The  effect  of  load  and  r on  stress  concen- 
tration factor.  c 
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figure  7.10,  required  an  average  running  time  of  10  min  each  on  the 
IBM/360/65  computer,  corresponding  to  $8.33  per  case,  with  N = 34 
nodes  per  case. 

7.6  Cone! usions 

A non-Hertzian  elastic  contact  problem  involving  an  unknown 
multiply-connected  contact  region  has  been  solved.  The  example 
problem  considered,  is  that  of  a pitted  sphere  in  contact  with  an 
unpitted  sphere.  The  axi symmetry  of  the  problem  enabled  us  to  use 
the  "simply-discretized  method"  with  a polar  coordinate  grid.  For 
problems  with  a lower  degree  of  symmetry,  it  had  been  found  in  earlier 
work,  that  a more  complicated  (and  less  efficient)  method  of  solution 
was  necessary  because  of  the  numerical  instability  of  the  equations 
generated.  It  may  be  appropriate  to  describe  the  equation  set  (7.24) 
as  "quasi-stable"  because  it  exhibits  dependence  on  the  location  of 
the  (N  + 1 ) th  field  point.  Through  experience  and  heuristic  reasoning, 
it  was  established  that  locating  the  additional  field  point  ( P ' in 
fig.  7.3)  at  the  outside  boundary  yields  a wel 1 -condi tioned  matrix. 

The  variation  of  the  SCF,  contact  region  and  peak  pressure 
Pmax  changes  in  the  pit  blending  radius  r*,  and  the  pit  edge 
radius  r*,  was  studied,  and  some  numerical  results  were  presented. 

The  numerical  solution  was  shown  to  converge  rapidly  with  a 
moderate  cell  density. 

The  principal  results  of  this  chapter  have  been  published 
by  the  International  Union  of  Theoretical  and  Applied  Mechanics  in 
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a joint  paper  by  Woodward,  Paul  and  Singh.  To  the  best  of  our 
knowledge,  this  is  the  first  published  solution  of  a multiply- 
connected  contact  region  problem  with  an  a priori  unknown 
contact  boundary. 


8.  CONCLUSIONS 


A general  method  of  solution  of  three  dimensional  frictionless 
conformal  contact  problems  has  been  presented.  Specifically,  two  con- 
formal examples  were  analyzed,  viz.,  the  case  of  an  elastic  sphere 
indenting  an  elastic  spherical  seat  and  the  case  of  an  elastic 
cylinder  indenting  an  elastic  cylindrical  seat.  The  necessary  influ- 
ence functions,  needed  for  solution  of  these  problems,  were  generated 
numerically  and  validated  with  analytic  solutions  wherever  possible. 

The  predicted  values  of  contact  stress,  load,  approach,  and 
contact  area  for  these  examples  is  in  close  agreement  with  Hertz's 
solutions  in  the  case  of  small  loads,  where  small  contact  regions 
occur.  For  larger  angles  of  contact,  the  load-approach  relationship 
obtained  for  the  sphere-seat  problem  was  found  to  compare  favorably 
with  the  experimental  results.  The  displacement  field  obtained  in  the 
analysis  of  an  elastic  sphere  in  contact  with  a conformal  rigid  seat 
was  found  to  be  reproduced,  with  satisfactory  accuracy,  by  a finite 
element  model  subjected  to  the  same  pressure  distribution. 

Therefore,  it  may  be  concluded  that  the  solution  obtained  by 
the  methods  of  this  dissertation  is  the  unique  solution  to  the  problem. 
Pressure  distributions  and  maximum  pressure  obtained  in  the  anlaysis  of 
conformal  cylinders  were  within  1.5  percent  of  the  values  predicted  by 
Persson  [1964].  The  values  of  load  vs.  contact  angle  also  agree  with 
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the  values  computed  by  Persson.  In  reviewing  the  above  results,  it 
is  concluded  that  the  method  presented  can  be  used  to  successfully 
analyze  three  dimensional  conformal  contact. 

In  addition  the  non-conformal  problem  involving  a pitted- 
sphere  contacting  a sphere  was  solved  for  a variety  of  pit  geometries. 
This  problem  is  of  mterest  because  it  has  a multiply  connected 
contact  region.  The  necessary  iteration  needed  to  converge  to  the 
true  boundaries  of  the  problem  was  established.  The  predicted  values 
of  contact  stress  were  essentially  Hertzian  away  from  the  pit  location 
and  the  stress  became  much  larger  in  the  vicinity  of  the  pit,  as 
expected;  the  stress  concentration  factor  was  found  as  a function  of 
pit  geometry  parameters. 

In  all  of  the  above  mentioned  axisymmetric  examples, 
axisymmetry  was  utilized  to  the  fullest  extent  in  the  discretization 
process  and  it  was  discovered  that  all  solutions  were  quasi-stable 
using  the  simply  discretized  method  of  Singh  and  Paul  [1974],  This 
has  not  been  observed  before. 

The  computer  costs  in  all  cases  were  minimal,  being  at  most 
$8. 33/for  the  34  field  point  model  for  the  pitted  sphere  examples. 

The  principle  conclusions  of  the  foregoing  results  may  be 
summarized  as  follows: 

1.  A general  method  of  solving  frictionless,  three  dimensional 
conformal  contact  problems  has  been  formulated. 
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2.  A method  by  which  numerical  influence  functions  may  be  generated 
is  presented  and  val idated  where  analytic  solutions  could  be 
obtained. 

3.  The  conformal  analysis  was  shown  to  be  in  close  agreement  with  the 
limiting  cases  of  Hertzian  contact  for  light  loads  and  with  other 
numerical,  analytic  and  experimental  analyses  of  conformal  contact 
problems. 

4.  A problem  involving  a multiply  connected  contact  region  was  solved. 

Future  research  should  be  directed  towards  applying  the  above 
method  to  the  solution  of  non-axi symmetric  problems  in  conformal  con- 
tact. Within  the  broad  area  of  elastic  contact  theory,  the  inclusion 
of  friction  and  dynamics  in  contact  theory  are  needed  areas  of 
investigation. 


APPENDIX  A 


DOMINANT  SINGULARITIES  IN  THE  STERNBERG  INFLUENCE 
FUNCTION  FOR  A POINT  LOAD  ON  A SPHERE 


Consider  a sphere  of  radius  R compressed  between  two  diamet- 
rically opposed  point  loads,  F.  The  displacements  u^  and  uQ,  as  shown 
in  figure  A.l,  are  derived  in  Lure  [1964]  and  are  given  by  equations 
(A.l)  and  (A. 2).  It  will  be  shown  in  this  appendix  that  the  dominant 
singular  terms  in  these  displacement  functions  are  those  of  the 
Boussinesq  influence  functions  for  a point  load  on  a plane. 

The  Sternberg  Influence  Functions  may  be  written  as  follows: 
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Fig.  A.l . Diametrically  Opposed  Point  Loads  on  a Sphere 


Fig.  A. 2. 


Coordinate  System  for  a Point  Load  on  a Plane 
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where  m is  the  reciprocal  of  Poisson's  ratio,  G is  the  modulus  of 
rigidity,  P2k  (cos  0)  are  the  Legendre  polynomials  in  cos  0 and  the  con 

stants  Ap  and  Bn  are  given  by 
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where 


a'  = n fn-/)  + (2  mi) 


m+ 1 


(A. 5) 


The  Boussinesq  displacement  functions  are  given  by 


Uy  = 


fp-n2) 
it  er 


(A. 6) 


and 
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where  v and  E are  Poisson's  ratio  and  Young's  modulus  respectively  and 
uz>  u^  and  r are  illustrated  in  figure  A. 2.  Consider  the  coordinate 
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system  in  figure  A.l  where  r = R0.  As  r tends  to  zero,  r will  be 
assumed  equal  to  distance  r along  the  straight  boundary  in  figure  A. 2. 
Also  0 in  equations  (A.l)  and  (A. 2)  will  be  replaced  by  r/R. 
Specifically,  it  will  be  shown  that 

j&vh  Ur  r ^2  ^A'8^ 

and 
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Now  consider  the  singular  terms  in  equation  (A.l).  They 
are  as  follows: 


F m-i l_ 

4rrGR  M % 


(A. 10) 


F_ 

4wGR 


(A. 11) 


4TTGR 


(■ 


m 


V . t e a \ 


It  can  be  seen  upon  inspection  that  all  other  terms  tend  toward  a 
finite  quantity.  and  definitely  tend  toward  » as  0 tends  toward 
zero  whereas  is  indeterminant  at  0 = 0.  It  will  be  shown  for  the 
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limit  of  r tending  to  zero  that 

S,  3 


(A. 13) 
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Expanding  S-j  , in  terms  of  0 we  find 
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As  0-*O  the  only  singular  term  is  the  first,  which  will  be  termed  S 


b 
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Now  consider  the  lim  of  the  ratio 

r+0  u 

m = v and  G = 2n^y  • 
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Recal 1 i ng  the  0 = r/R, 
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Therefore  the  singular  term  S-j  corresponds  directly  to  u . 
Investigating  the  limit  of  the  ratio  S-j/S^,  where 
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S,  _ m(m-i)  <4*  f 
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it  can  be  shown,  using  1'Hopital's  rule  that 
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Thus,  even  though  the  singularity  is  present  in  the  displacement 
ur>  it  is  weak  compared  to  the  singularity  S-j . 

Finally,  examining  S^,  using  1'Hopital's  rule  it  can  be 
shown  that 


(A. 20) 


In  other  words,  is  not  a singular  term,  but  tends  toward  zero  at 
the  load. 

From  the  above  analysis  it  can  be  seen  that  the  dominant 
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singularity  in  u^  is  and  furthermore,  in  the  limit  as  r tends  toward 

zero  S-j  = uz-  Therefore,  it  may  be  concluded  that  in  the  same  limit 

u = u . This  is  illustrated  in  table  A.l.  As  0 decreases,  the  per- 
r z ' 

cent  difference  between  u and  u also  decreases. 

r z 

Now  consider  the  displacements  u()  and  ufc.  The  three  terms  in 
equation  (A. 2)  which  are  singular  or  indeterminant  are 
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where  k2  = 1H|M  . 

T-j  is  singular  while  T2  and  T^  are  indeterminant.  The  other  term  in 
equation  (A. 2)  involving  dPo^cos  0)/d6  is  known  to  be  zero  at  0 equal 

Examining  the  ratio  T^/u^  it  can  be  shown  that  for  r =R0 


to  zero. 


196 


TABLE  A. 1 

COMPARISON  OF  BOUSSINESQ  INFLUENCE  FUNCTION  TO  THE 
INFLUENCE  FUNCTION  FOR  A POINT  LOAD  ON  A SPHERE 


0[Deg.] 

u [in.]a 
rL  J 

uz[in.] 

7o  Difference 

10.00 

1 .44963 

1 .65964 

12.654 

1 .00 

16.4732 

16.5964 

0.74222 

0.10 

165.917 

165.964 

0.02806 

0.01 

1659.65 

1659.64 

-0.00179 

R = 10",  F = 30  x IQ7  lb,  E = 30  x 106  psi , v = 0.30 


a51  Terms  taken  in  Legendre  Polynomial 


_7 ; 

ut 
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Usi ng  1 'Hopital ' s rule. 
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Therefore  in  the  limit  as  0,  or  r,  tends  to  zero  = T-j . 

Consider  term  in  the  limit  as  0 goes  to  zero,  and  using 
1 'Hopital ' s rule 
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This  proves  that  although  T^  was  indeterminant  it  is  not  a singular 
term  in  the  limit  as  0 tends  to  zero. 

Finally,  consider  T^.  T^  may  be  rewritten  in  the  form 
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which  is  indeterminant  as  8 tends  to  zero.  Using  l'Hopital's  rule 
and  reducing  the  results 


jLlfflr 


T O'  " ***  i 

S \-+0  8rrG  R no1  \/U*6  /*if(l+A*i) 


Aiy?&  ‘kntrz 

ao*  e(n-z*i£) 


(A. 28) 


Therefore  is  not  a singularity. 

From  the  above  analysis  it  was  shown  that  the  only  singularity 
in  the  function  u^  is  T-j  and  that  this  is  equal  to  u^  as  0 tends  to 


zero . 


APPENDIX  B 


DISPLACEMENTS  ON  A CYLINDER  UNDER  TWO 
DIAMETRICALLY  OPPOSED  LINE  LOADS 

Consider  a cylinder  under  two  lines  loads  F as  shown  in 
figure  B.l.  Muskhel i shvi 1 i [1963]  has  shown  that  for  the  case  of 
plane  strain  the  displacements  in  the  x and  y directions  respectively 
are. 
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Fig.  B.2.  Geometry  for  Diametrically  Opposed  Line  Loads  on  a Cylinder 
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where  A and  p are  the  Lame  constants. 

Now  consider  the  case  where  the  two  line  loads  are  diametri- 
cally opposed  (figure  B.2).  The  following  relations  may  be  derived 
for  the  geometry  illustrated  in 
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Substituting  in  the  relations  of  (B.2)  into  equations  (A.l)  and  (A. 2), 

the  displacements  u and  u become, 

K x y 

Ux  T F K,  [JIm  fa*.  I f I + 604  & | (b.4) 
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where 

* CA+v)  _ z(i-v') 

a.  ~ r — — — ( b .6 ) 

2 TT/U.  (h/U.)  TTE 

I (i-n>)(i-zv) 

K>  - — , ~ r - r~zr — — (b.7) 

4fA  +m)  2 f 

E and  v are  Young's  modulus  and  Poisson's  ratio  respectively.  In 

polar  coordinates  the  radial  and  tangential  displacements,  notated  u^ 

and  u.  respecti vely,  computed  as  functions  of  u and  u are 
t x y 

Ufr  • — Ux  Q — My  0 

(Jq  s — AM?*  & + (Jy  Q 

Substituting  equations  (B.4)  and  (B.5)  into  equations  (B.8)  and  (B.9), 

u and  uQ  become 
r 0 


(B.8) 
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and 
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APPENDIX  C 


SINGULARITIES  IN  THE  INFLUENCE  FUNCTIONS  FOR  A 
CYLINDER  UNDER  TWO  DIAMETRICALLY 
OPPOSED  LINE  LOADS 

Consider  a cylinder,  of  radius  R,  under  two  diametrical ly 
opposed  line  loads  F,  as  shown  in  figure  C.l.  The  displacements  u^ 

and  u.  as  derived  in  appendix  B are 

0 

(Jr  - - K,F  JUi  -km  J + l]  + jt  Fa**  6 (c.i) 
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Fig.  C.l.  Two  Diametrically  Opposed  Line  Loads  on  a Cylinder 


Fig.  C.2.  Line  Load  on  an  Elastic  Half  Space 
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E is  Young's  modulus  and  v is  Poisson's  ratio.  Also  consider  the 
Flamant  solution  for  a line  load  acting  on  a half  space  as  illustrated 
in  figure  C.2.  As  presented  in  section  4.3,  the  displacements  uz  and 
ut  are 


tte  * r 


(\+v)F 

THE 


(C.5) 


and 


(l  + v)(l-2V)  F 

ZE 


(C.6) 


it  will  be  shown  that  the  singularity  in  the  displacement  u^  is  the 

same  as  that  of  u near  the  load.  Furthermore,  the  limit  of  uA , as  0 

z 0 

tends  to  zero,  will  be  shown  to  be  u^. 

As  illustrated  in  figure  C.2,  r represents  the  distance 

between  the  applied  load  and  the  point  Q where  the  displacements  uz 

and  u^  are  calculated.  Let  us  define  a coordinate  r on  figure  C.l 

such  that  r = R0.  Replacing  0 with  ^ in  equation  (C.l),  it  will  be 

shown  that  u - u as  r tends  toward  zero.  Likewise,  uA  will  be  shown 
r z 0 

to  approach  u^  as  r tends  to  zero. 

In  terms  of  the  coordinate  r,  the  ratio  u^/uz  may  be  written 


as 
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where 
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Taking  the  limit  as  r tends  to  zero  of  equation  (C.7)  and  using 
1 ' Hopi tal 1 s rule. 
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The  second  term  in  (C.9)  clearly  goes  to  zero  for  r equal  to  zero 
while  the  first  and  third  terms  are  indeterminant. 

Expanding  cot  in  the  first  term  of  equation  (C.9),  and 
taking  the  limit  as  r tends  to  zero 
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~ j (C.10) 

Taking  the  limit  of  the  third  term  in  equation  (C.9)  and  using 
1 'Hopital 's  rule. 


r ..  i r , • * r 
ZR  A£c-  ZR  ***  R 
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Hence,  the  only  nonzero  term  in  equation  (C.9)  in  the  limit  as  r tends 
to  zero  is  the  first  and  it  tends  toward  1.  Therefore  in  the  limit 
as  r tends  to  zero,  ur  = u^. 

Now  consider  the  limit  of  un/u.  . Both  uQ  and  u.  are  finite  at 

U L u L 

r equal  to  zero.  In  equation  (C.2),  the  first  term  can  be  shown  to  go 
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to  zero  for  small  r.  uQ,  therefore,  tends  to  the  value  -K^  for  small  r. 
This  is  exactly  the  value  of  u^  for  all  r.  Thus, 


Jb*.  Jdl.  , / (C.12) 

r—o  Ut  1 


and  uQ  = u.  when 


r tends  to  zero. 


APPENDIX  D 


DERIVATION  OF  SURFACE  DISPLACEMENTS  FOR  A 
CYLINDRICAL  CAVITY  UNDER  TWO 
DIAMETRICALLY  OPPOSED 
LINE  LOADS 

Timoshenko  and  Goodier  [1970]  derived  the  stress  functions 
^ and  (|>  corresponding  to  an  elliptic  hole  with  uniform  pressure  p 
applied  on  two  diametrically  opposed  segments.  Consider  a unit  circle 
in  the  c plane  and  the  mapping  function, 

£•  co(s)  <D1> 

where 

<»(r)  = R(  T+  f)  (d.2) 

a )(c)  maps  the  unit  circle  in  the  £ plane  into  an  ellipse  in  the  Z 
plane  with  semi  axes 


CL  r R(lfm) 

(D.3) 

b - R(!-m ) 

For  m = 0,  the  mapping  becomes  a circle  of  radius  R.  Now  consider  an 
elliptic  hole  with  a pressure  distribution  applied  as  shown  in 
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figure  D-l . Timoshenko  has  shown  that  the  stress  functions  in  the  c 
plane,  corresponding  to  the  above  loading,  are 


and 


+ 


(D.4) 


where  a-| , and  correspond  to  the  mapped  points  Z-j  and  in  the  Z 
plane.  Letting  m = 0 for  the  case  of  a circular  hole,  the  mapping 
function  becomes 


? oj (■?)  = /?  ? 


V i'i 

Fig.  D.l.  Elliptic  Hole  with  Internal  Pressure 


ie 
U , 


Fig.  D.2.  Diametrically  Opposed  Line  Leads  on  a Cylindrical  Cavity 
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and  equations  (D.4)  and  (D.5)  may  be  written  in  the  Z plane  as 
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The  displacements  u^  and  uQ  are  related  to  the  stress  functions  by 


2/p(ur+iu>)=  et6  <t>'(i)-yn)  J 


(D.8) 


where  p is  the  modulus  of  rigidity,  v is  Poisson's  ratio  and 


(D.9) 
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for  plane  stress  or 

X - J-4j/  (DJ0 

for  plane  strain.  (See  figure  D.2.) 

In  order  to  find  the  correct  displacements  for  concentrated, 
diametrically  opposed  loads,  F,  let  p vary  as  F/|Z^  - | and  take 

the  limit  as  |Z^  - |+0.  Hence  for  concentrated  loads  the  stress 

functions  become 
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The  quantities  Z,  and  Z^  may  be  written  as 

z = Rei@ 
z,  --  S e ** 

?,  = Re'1* 


(D.l  3) 


Substituting  equations  (D.13)  into  equations  (D.ll)  and  (D.12),  the 
functions  ip  and  <f>  may  be  written  in  terms  of  0 and  a.  The  limit  of 
|Z-j  - Z.j  | is  now  the  limit  as  a->0.  Taking  that  limit,  the  functions 
ip  and  (j)  become 
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where  Z = Zn  = R. 
o U 

The  displacements  u^  and  uQ  may  be  found  by  substituting 
equations  (D.14)  and  (D.15)  into  equation  (D.8).  Separating  out  the 
real  and  imaginary  parts,  u and  ufi  are  found  to  be 
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For  plane  strain  these  displacements  may  be  written  in  terms  of  E and 


v as 
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APPENDIX  E 


SINGULARITIES  IN  THE  INFLUENCE  FUNCTIONS  FOR  A 
CYLINDRICAL  CAVITY  UNDER  TWO  DIAMETRICALLY 
OPPOSED  LINE  LOADS 


Consider  a cylindrical  cavity  of  radius  R,  under  two  diamet- 
rically opposed  line  loads  F,  as  shown  in  figure  E.l.  u^  and  uQ  as 
derived  in  appendix  D are 
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E is  Young's  modulus  and  v is  Poisson's  ratio.  Also  consider  the 
Flamant  Solution  for  a line  load  acting  on  a half  space  as  illustrated 
in  figure  E.2. 
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Fig.  E.l.  Two  Diametrically  Opposed  Line  Loads  on  a Cylindrical  Cavity 


Fig.  E.2.  Line  Load  on  an  Elastic  Half  Space 
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As  presented  in  section  4.3,  the  displacements  uz  and  u^  are 


ZF(\-VX)  A o!  F(t+») 


(E.5) 


and 


(\+v)(i-iv)  F 
_____ 


(E.6) 


It  will  be  shown  that  the  singularity  in  the  displacement  is 

the  same  as  that  of  u near  the  load.  Furthermore,  the  limit  of  uQ, 

as  6 tends  to  zero,  will  be  shown  to  be  u^. 

As  illustrated  in  figure  E.2,  r represents  the  distance 

between  the  applied  load  and  the  point  Q where  the  displacement  uz  and 

ut  are  calculated.  Let  us  define  a coordinate  r on  figure  E.l  such 

that  r = R0.  Replacing  0 with  ^ in  equation  (E.l),  it  will  be  shown 

that  u = u as  r tends  to  zero.  Likewise,  u„  will  be  shown  to 
r z 0 

approach  u^  as  r tends  to  zero. 


In  terms  of  the  coordinate  r,  the  ratio  ur/uz  may  be  written 


as 


Ui  "A  ] 


(E.7) 


where 
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L . - <11^1 

K*  ' ire 


(E.8) 


Taking  the  limit  as  r tends  to  zero,  of  equation  (E.7)  and  using 
1 * Hopi tal 1 s rule 


£ Kx 
R 


( E.  9) 


The  third  term  in  equation  (E.9)  clearly  goes  to  zero  as  r tends  to 
zero,  while  the  first  and  second  terms  are  indeterminant. 

Expanding  cot  ^ in  the  first  term  in  equation  (E.9)  and 
taking  the  limit  as  r tends  to  zero. 


jur. 

R 


- — /£  _£ 
t-^o  R \ r 3R  4-SR3 


(E.10) 


Taking  the  limit  of  the  second  term  in  equation  (E.9)  and  using 
1 'Hopitals  rule 
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= O (E.ll) 


Hence  the  only  nonzero  term  in  equation  (E.9)  in  the  limit  as  r tends 
to  zero  is  the  first  and  it  tends  to  1.  Therefore  in  the  limit  as 


r tends  to  zero,  ur  = uz* 

Now  consider  the  limit  of  Ug/u^.  Both  uQ  and  are  finite 

at  r equal  to  zero.  In  the  first  term  of  equation  (E.2)  can  be  shown 
to  go  to  zero  as  r approaches  zero.  uQ,  therefore,  tends  to  the  value 
-K^Eforsmall  r.  This  is  exactly  the  value  of  u^  for  all  r.  Thus, 


\r-+0  Uf 


Am  = I 


(E.12) 


and  u0  = u».  when  r equals  0. 


APPENDIX  F 


DERIVATION  OF  THE  PROFILE  FUNCTION  FOR  CONFORMAL 


CONTACT  OF  A SPHERE  AND  SPHERICAL  SEAT 


Consider  a sphere  of  radius  R^  in  contact  with  a spherical 
seat  of  radius  R^  at  a point  0 as  illustrated  in  figure  F.l.  It  is 
assumed  that  point  A,  located  at  ip  on  the  sphere,  will  contact  point 
B on  the  seat,  located  at  <j>.  The  distance  between  A and  B,  denoted 
f,  is  the  profile  function  for  these  points.  The  value  of  the  pro- 
file function  in  terms  of  <J>  and  ip  is  derived  below. 


Vectors  r-j  and  r^  are  defined  such  that  r-j  extends  from  0 to 
A and  r2  extends  from  0 to  B.  The  quantities  w-|  and  define  the 
angles  between  the  x axis  and  r-j  and  r^  respectively.  From  geometry 
of  isosceles  triangles  it  can  be  shown  that 


'V 


2 


(F.l) 


and 


4> 


2 


(F.2) 


Furthermore,  from  geometry 


/ r,  I - Z Avn  y 


(F.3) 
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Fig.  F.l.  Profile  function  for  a sphere  and  spherical  seat 


A 


Fig.  F.2.  Profile  function,  f,  relative  to  mean  radial  direction,  r 


in. i --  zrz 
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(F.4) 


The  components  in  the  x and  z directions  of  r-|  amd  r?  may  be 


expressed  as 


r,%=  fi,  Ain't 
If  l ] g * 2 Rt  2 ^r* 

(F.5) 

^2%  ~ &z  $ 
fzi  -2^2.  J" 

The  values  of  r,  , r,  , r0  and  r0  are  merely  the  cartesian 
lx  lz*  2x  2z  J 

coordinates  of  points  A and  B.  The  value  of  | f | is  therefore  defined 
by  the  distance  between  these  two  points  or 

I fl=  -s/f  n%  - 'i-x  )l  + tt*-***)* 

where  r,  , r«  , r,  and  r0  are  given  by  relations  (F.5)  in  terms  of 
lx  2x  lz  2z  3 J 

ip  and  4). 

Now  consider  the  mean  radial  direction  which  forms  an  acute 

<f>  + <J>  , 

angle  of — ~ — "with  the  z direction  as  shown  in  figure  F.2.  The  angle 

between  AB  and  the  mean  radial  direction  is  £ while  the  acute  angle 
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between  AB  and  the  z axis  is  labeled  y such  that  £ + y =(4)  + 40/2. 
The  value  of  £ may  be'  determined  using  the  following  relations 


(F.7) 


(F.8) 


APPENDIX  G 


INTEGRATION  OF  THE  BOUSSINESQ  INFLUENCE  FUNCTION 
OVER  AN  ANNULAR  ELEMENT 


It  is  required  to  evaluate  the  integral 


zav  -Jfj 


(x-x'JVs- ay 


(G.l) 


over  an  annular  element  where  (x,y)  represents  the  location  of  a field 
point  C along  the  center  line  of  the  annular  segment  as  illustrated  in 
figure  G.l.  The  annular  element  has  inner  radius  r^ , outer  radius  r^ 
and  the  field  point  is  located  as  radius  r . The  sides  of  the  element 
are  defined  by  angle  9 measured  from  the  center  line. 

Now  consider  the  right  triangle  as  shown  in  figure  G.2.  The 
integral  in  equation  (G.l)  over  the  right  triangle  is  given  in  Lure 
[1964]  as 

X ~ J t&M,  ( + y ) (G.2) 

This  result  will  be  used  to  approximate  the  integral  in  equation  (G.l) 
over  the  annular  region  in  figure  G.l. 

Let  the  annular  region  of  figure  G.l  be  divided  into  six  right 
triangles  as  shown  in  figure  G.3.  They  are  as  follows. 
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* 


Fig.  G.l.  Annular  Segment  with  Field  Point  along  Center  Line 


Fig.  G.2.  Right  Triangle  with  Field  Point  as  a Vertex 
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Fig.  G.3.  Annular  Region  Subdivided  into  Six  Right  Triangles 
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The  calculation  of  the  integral  (G.l)  over  each  triangle  may  be  per- 
formed using  relation  (G.2).  The  integrals  over  each  of  the  triangles 

defined  by  relations  (G.3)  will  be  termed  I-j,  1^  corresponding 

to  the  integrals  over  triangles  1,  2 6 respectively.  Thus 

li'JiUtto+d*  |‘)  (G4) 

where  <J>.  i = 1,6  are  illustrated  in  figure  G.3.  The  values  of  <j>.  are 
all  functions  of  r-| , r^,  rc  and  0 defined  below. 


a - I 


(G.5a) 
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c (rc-  r,  )/t^i  f i 

| -===  = } 

W 'i'W©  ♦O'icoie-r;)1  j 


(G.5b) 


03  = mtJ&m,  ^ 


(rc  cae-r,) 

(*c  & ) 


■) 


4 t (£.4^)  J 


4 - JC  „ 


(G.5c) 


(G.5d) 


(G.5e) 


^ ~ rr  - 0a  - -a4  - 

The  values  of  JL  can  then  be  defined  by 

1,  - Jli^  CH  - (rc-n  )a>d  f 

/.  -J4  ~ c£  - rc  a**.  © 


(G.5f) 


(G.6a) 


(G .6b) 


Xs  - A - s^fc-4**0),+(V£C<*0Jt 


(G.6c) 


The  total  integral  ly  over  the  annular  element  can  be  approximated 
by 


XT  * ^ {i/<  A Ah  f f * # > / A4.  rf ♦ £ )}  <G'7) 


APPENDIX  H 


DERIVATION  OF  A CONTACT  CRITERION  FOR 
CLOSELY  CONFORMING  SPHERES 

Consider  a sphere  of  radius  R^  interpenetrating  the  surface 
of  a spherical  seat  of  radius  R^  by  an  amount  6.  For  the  problem  at 
hand  6 corresponds  to  the  rigid  body  approach  due  to  some  unknown 
applied  force  F,  on  the  sphere.  Figure  H.l  illustrates  the  sphere 
and  seat  and  pertinent  notations  for  the  discussion  to  follow. 

Point  A is  located  on  the  surface  of  the  sphere  at  coordinate 
ip.  Point  B is  located  on  the  seat  between  point  A on  the  sphere  and 
02*  the  center  of  the  seat.  The  radial  gap  between  A and  B,  denoted 
by  f,  is  equal  to  the  sum  of  the  elastic  radial  displacements  w-j  and 
W2*  on  the  sphere  and  seat  respectively,  in  the  deformed  state.  From 
the  geometry  of  the  problem  the  following  relationships  may  be  stated: 

CL  A = 0}C  + CA  - (e+<5)/6w  d +/?,  tot  (V-$) 

(H.l) 

< \&~Ri  (H.2) 

6 A » Oz/A  -OtS  u(e-tl»)£#A<fl+R,yt4ri('f-$)'Rz 

(H.3) 

where  e = R^  - R-j  • 
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Fig.  H.l.  Sphere  interpenetrating  a spherical  seat 
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Now  consider  the  case  where  ~ R-j  and  6 « R-| , hence  ip  ~ <t>. 
With  these  assumptions  equation  (H.3)  may  be  written  as 


BA  ^ (H . 4) 

Recalling  that  BA  = w-j  + 

l V,+^2  t 6 (1-jCsn  f)  y - O (H.5) 

Equation  (H.5)  is  an  approximate  contact  criterion  based  on  radial 
displacements. 


APPENDIX  I 


RELATIONSHIP  BETWEEN  THE  ELASTIC  CONSTANTS 
IN  PLANE  STRESS  AND  PLANE  STRAIN 

Given  the  displacement  and  stress  fields  in  a state  of  plane 
stress  in  terms  of  the  elastic  constants  E and  v,  the  equivalent 
fields  for  the  identical  problem  in  plane  strain  may  be  found  by 
substitution  of  E for  E and  v for  v where 


and 


(1.1) 


(1.2) 


This  can  be  verified  by  substitution  of  (I  • 1)  and  (1*2) 
into  the  stress-strain  relations  for  plane  stress.  The  stress-strain 
relationships  for  a linear  isotropic  material  are 


£% s £ i ) 

- Y «r% -»<rx) 

v . T 

" £ l** 


(1.3) 
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Substituting  equations  (1*1)  and  (I  • 2)  into  (I  • 3)  yields 


Equations  (I  • 4)  are  the  stress-strain  relations  for  plane 

strain. 

Similarly,  a solution  in  plane  stress  may  be  obtained  from 
a solution  in  plane  strain  by  substituting  E'  for  E and  v'  for  v 
where 


<fx  s y|W)o;-  j 

- jiO'Vl)c r,-VO+v)0-K  J 


(1.4) 


' 


(l+v)z 


(1.5) 


(1.6) 


Substituting  relations  (I  • 5)  and  (I  • 6)  into  equations  (I  • 4)  one 
obtains  equations  (1*3)  for  plane  stress. 


APPENDIX  J 


HERTZIAN  FORMULAS  FOR  A SPHERE 
INDENTING  A SPHERICAL  SEAT 


For  the  case  of  a sphere  of  radius  R^  indenting  a spherical 
seat  of  radius  R2»  Hertz's  theory  predicts  the  following  relation- 
ships between  load  F,  approach  6,  and  the  radius  of  the  contact 
region  a (see  Timoshenko  and  Goodier,  1970  pp.  409-14): 

I 


a = 


£ - 


3rr  k Pl  F 

F*  k*<Rt- B,) 


(j.i) 


16 


r,r2 


(J.2) 


where 


, i -u(*  i-vS 

R - — + 


(J.3) 


W Ft  ITFj 

v-j  and  \>2  are  Poisson's  ratio  for  the  sphere  and  seat  respectively 

and  E^  and  E^  are  the  respective  values  of  Young's  modulus. 

The  pressure  at  a radius  of  r from  the  center  of  the  contact 
region  is  given  by 


P(r )* 


3 F 

2VQ.1 


I- 


(J.4) 
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These  formulae  are  only  valid  for  contact  regions  which  have  dimen- 
sions small  compared  to  both  R-j  and  R^. 


APPENDIX  K 


P 


HERTZIAN  FORMULAS  FOR  A CYLINDER 
INDENTING  A CYLINDRICAL  SEAT 


For  the  case  of  a cylinder  of  radius  R-j  indenting  a cylin- 
drical seat  of  radius  R^,  Hertz's  theory  predicts  the  following 
relationships  between  load  F per  unit  length  and  the  half  width  of  the 
contact  region  b (see  Timoshenko  and  Goodier,  1970  pp.  418-19): 

a 

( K.  1 ) 


b = 


4Fk  R,Rl 

Rt-R, 


bl  (Pi-R.) 
4 k P,  Pt 


(K.2) 


where 

k = 

T£ ; TTft 


(K.3) 


and  are  the  Poisson's  ratio  of  the  cylinder  and  seat  respec- 
tively while  E^  and  E^  are  the  respective  values  of  Young's  modulus. 

The  pressure  distribution  at  a distance  of  r from  the  center 
line  of  contact  region  is  given  by 
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These  formulae  are  only  valid  for  b « R-j . The  approach  for  cyl- 
inders in  contact  or  the  contact  of  a cylindrical  seat  is  predicted 
to  be  infinite  by  Hertzian  theory  which  is  clearly  not  possible.  For 
an  explanation  of  this  inadequacy  and  the  derivation  of  an  appro- 
priate formula  for  cylinlers  in  line  contact  see  Singh  [August  1974]. 


APPENDIX  L 


DERIVATION  OF  A CONTACT  CRITERION  FOR  CONTACT  OF  A SPHERE 
AND  SPHERICAL  SEAT  WITH  A CONSTRAINED 
DISPLACEMENT  FIELD 

Consider  an  elastic  sphere  of  radius  R-j  indenting  an  elastic 

spherical  seat  of  radius  R^.  For  a given  force  applied  to  the  sphere, 

the  bodies  will  approach  by  6 and  the  contact  area  will  extend  to  \p 

rr  max 

on  the  sphere.  It  will  be  assumed  that  a point  located  at  ip  on  the 
sphere  will  contact  a point  at  4>  on  the  seat  and 

Ri  _ Aton.  0 
Ri  a*a* 

This  is  physically  equivalent  to  requiring  points  with  the  same  x 
coordinates  as  shown  in  figure  L.l  to  contact  after  deformation.  Only 
displacements  in  the  z direction  will  be  considered  in  the  criterion. 
Consider  the  sphere  interpenetrating  the  seat  (although  this  is 
physically  impossible)  by  an  amount  6.  The  distance  between  points  A 
and  B is  labeled  f. 

Noting  the  geometry  at  hand  and  the  constraint  equation  (L.l) 
the  gap  f may  be  written  as 

■f  - (ftx-R,)4  & x/n  'f  -Ri  4>  (l.2) 
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Fig.  L.l.  Sphere  interpenetrating  a spherical  seat 


interpenetrate^ 


w,  AND  ARE  DiML  N^lONl  ESS 
Fig.  L.2.  Detail  of  interpenetration  surfaces  of  sphere  and  seat 
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This  gap  must  be  closed  by  the  dimensionless  elastic  displacements 
w1  and  v) ^ on  the  sphere  and  seat  respectively,  (see  fig.  L.2)  Thus 

R,  K tRx  ^2  = (&i  ~R\)  t & ~ (L.3) 


Writing  equation  (L.3)  for  points  on  the  outer  boundary 

J^iWi  +R%  Vvjj  Ri-Rt  + 4 ~(Ri £0* ^ '%*#%) 

(L.4) 

Equation  (L.4)  may  be  solved  for  6.  Substituting  the  resulting  value 
of  6 into  equation  (L.3)  and  rearranging  terms  gives 

( tv,  -iv( a)R,  +(wz  -w°)  r%  - 

~ /?>  (L.5) 


Writing  equation  (L.l)  for  points  on  the  outer  boundary 

R\  Am 

— - — - 1 - — (L.  6) 

R i Am 

Dividing  equation  (L.5)  by  R-j  and  substituting  in  relation  (L.6)  the 
final  expression  for  the  contact  criterion  becomes 


Am.  t/v\a* 
Am  ftrrtAA 


( Wrwf)  + (iv*-w/) 

~ ZZfa/4*4?  Am*  )*  (** 


(L.7) 
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The  above  expression  (L.7)  is  the  contact  criterion  used  in  the 
analysis  of  Goodman  and  Keer  [1965]  who  omit  the  bar  over  the  dimen- 
sionless quantities  , w^. 


APPENDIX  M 


DERIVATION  OF  CONSTRAINED  DISPLACEMENT  FIELD  FOR 
SURFACE  POINTS  OR  AN  ELASTIC  SPHERE 
CONTACTING  A RIGID  SEAT 


Consider  an  elastic  sphere  of  radius  R^  in  contact  with  a 
rigid  spherical  seat  of  radius  R^.  It  is  known  a priori  that  the 
contact  surface  will  be  of  radius  R^.  Assume  that  all  points  on  the 
surface  of  the  sphere  move  in  a direction  parallel  to  the  line  of  the 
applied  load.  For  a given  contact  angle  ^max>  the  displacement  field 
is  then  uniquely  determined  and  will  be  derived  below. 

The  displacement  of  surface  points  on  the  sphere  consists  of 
a rigid  body  translation  6 and  an  elastic  displacement  field  uz  (ip) . 
Consider  a sphere  which  has  undergone  the  rigid  body  translation  so 
that  its  surface  interpenetrates  the  seat  as  illustrated  in  figure  M.l. 
The  elastic  displacement  of  point  A on  the  sphere  is  then  defined  by 
uz  in  accordance  with  the  given  assumption.  Taking  vector  components 


_A 


A 


of  OB,  O'A  and  00'  the  magnitude  of  uz  can  be  derived  to  be 


|ua|  ~ - loeljUH  Y i"  /oo'l  + I o' A 1 4&L  $ (Ml) 


244 


Fig.  M.l.  Elastic  sphere  interpenetrating  a rigid  seat 


i 


Fig.  M.2.  Components  of  surface  displacement  on  sphere 
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which  can  be  expressed  as 

lUtl  - Rz  (l-sCt±  <p)  - R,  fl  -44*  Y)  ■*&  (H.2) 

Let  a r - t coordinate  system  be  defined  at  point  A such 
that  ? is  directed  radially  inward  on  the  sphere  and  t is  directed  tt/2 
clockwise  of  ?,  Then  the  components  of  tTz , w^r  and  u-jt,  may  be  com- 
puted by 


w,*  )uzl  M+V 


(M.3) 


and 


u;  * / Uj  / aut 


(M.  4) 


The  initial  assumption  requires  A and  B to  be  located  such  that 


thus 


Rt  Aiw  7^  " Rt  aJm. 


(M.  5) 


(M.6) 
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Combining  equations  (M.3),  (M -4) , and  ( M . 6 ) one  finds 


W,  - Rt  ^1  -x&i  ( -j^r  A^y  )J  j ao<  y 

~ R,  ( I - A64  Y)  AAH  T + 6 604  Y 


(M-7) 


and 


W/ r i^/  ~ta*t 


(M.  8) 


du 


The  derivative  may  be  expressed  as 


- A*  AO<  £ AAAAon  f Ami  y)  j A&4.  Y 


- R,AiM*Y+Ri  co*Y  - R,  (i-am '})/**'}  +S oh.Y 


This  latter  derivative  is  useful  in  computing  the  surface  strain 
£#' 


APPENDIX  N 


COMPUTATION  OF  THE  ORIENTATION  OF  TANGENTIAL  COMPONENTS 
OF  DISPLACEMENT  ON  SPHERICAL  SURFACES 

Consider  the  spherical  surface  of  radius  R illustrated  in 
figure  N.l.  The  surface  is  subjected  to  a load  F at  point  B located 
at  (ipg , 6^).  The  tangential  displacement  ut»  at  a given  point  A 
located  at  (iK , 0),  lies  along  the  tangent  to  the  great  circle 
connecting  A and  B and  is  directed  away  from  B for  positive  values. 

The  angle  between  the  positive  ut  direction  and  the  tangent  at  A which 
lies  in  the  n-|  - y-j  plane  is  defined  at  t.  4>  is  the  angle  measured 
between  OA  and  OB.  For  purposes  of  integration  of  the  tangential 
displacement  influence  function  it  is  desired  to  find  x in  terms  of 

^A»  ^0  ^0- 

Consider  the  portion  of  the  spherical  surface  A M B as 
illustrated  in  figure  N.2.  The  spherical  angle  at  M is  3g  and  the 
opposite  angle  is  4>.  The  adjacent  angles  to  M are  ^ and  4^.  The 
spherical  angle  at  A is  x.  From  the  law  of  cosines  for  spherical 
trigonometry 


/CO40  - iC&iYg,  MHTa  *MMYe  (n-1) 
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A 


Fig.  N. 1.  Orientation  of  tangential  displacements  due  to  a 
point  load  on  a sphere. 


M 


B 


Fig.  N.2..  Spherical  triangle 
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The  law  of  sines  for  spherical  trigonometry  states: 


.W  T = 

r — (N.2) 

6 

combining  equations  (N.l)  and  (N.2)  and  taking  the  inverse  of  sin  t, 


T - OAX44m 

fii  Aim  Yb 

^cn  Act  V8  otn  % *Au«  Yg  Aim 

(N.3) 

The  above  result  can  also  be  derived  with  cartesian  vectors. 


APPENDIX  0 


INTEGRATION  OF  INFLUENCE  FUNCTION 
FOR  CYLINDRICAL  GEOMETRIES 

The  influence  functions  for  a cylinder  under  two  diametrically 
opposed  line  loads,  normal  to  the  surface  as  derived  in  section  4.7 
are 

<r,<  Y,ri  «..,Os  - kI  lM.fa»Cr~  + /}  + 

if  > 

f —2-  AtM.  ( V - y) 


(0.1) 


- !*Lco4(r'-Y) 

7 T 


(0.2) 


where  G-j  (if> , ip'  , E-j , ) is  the  influence  function  for  the  displace- 

l I 

ment  of  a point  at  ip  due  to  a load  at  ^ (^  - ip)  and 

H.|  ip' , E.j , v-|)  is  the  influence  function  for  the  tangential 
displacement-  E^  and  are  elastic  constants  of  the  cylinder  and  the 
constants  and  K l,  are  given  by 


K,  * 


Ki-- 


20-*i) 
7 TEi 

l£Ti 


(0.3) 

(C-Ijl) 

(0.4) 
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Similarly  for  the  cylindrical  seat  under  two  diametrically  opposed 
normal  line  loads  the  influence  functions  as  derived  in  section  4.8 


are 


s £*  /am.  - 

*’*  K\  Aim  ( 


(0.5) 


and 


Ht (t t: Vt)Et)  ‘ (+'- *) U 7 

- /<l  4*1  cv-t) 


(0.6) 


where  k^2  and  are  given  by  equations  (0.3,  0.4) 

G£  represents  the  radial  displacement  influence  function  and 
is  that  for  the  tangential  displacement  on  the  cylindrical  seat. 
Consider  the  loading  p (ip 1 ) on  the  cylinder  between  4>,  and 
\l>2  (4*2  > ) as  shown  in  figure  0.1.  p ( 4>  ' ) is  a constant  pressure 

P and  \p  is  located  such  that  ip  < 4^  (see  fig.  0.1).  The  displacement 
w-j  at  \p  due  to  p ( 47 ' ) can  be  computed  as 


The  tangential  displacement  un  due  to  load  p ( ip  ' ) can  be  expressed  as 


(0.7) 


(0.8) 
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Fig.  0.1.  Pressure  loading  on  a cylinder 


Fig.  0.2.  Pressure  loading  on  a cylindrical  seat 
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Similarly  on  the  cylindrical  seat  the  radial  and  tangential  displace- 
ments at  p due  to  p (<f) ')  are 


r* 

w2:  PR,  J 6,(ApyiiEt)d^' 


(0.9) 


and 


Ut  = PR i J 


(0.10) 


respectively  (see  fig.  0.2) 

The  integrals  in  equations  (0.7),  (0.8),  (0.9)  and  (0.10) 
may  be  evaluated  analytically.  In  each  case  the  first  term  in  G-j , 
H.| , and  can  be  integrated  by  parts  while  the  integration  of 
the  second  terms  is  trivial.  The  indicated  integration  results  in 

w,  - pr,  ^ -k,an(v'-r)f^x 


(0.11) 


-jinAtM  (Y-Y)) 

, lY  = Yt 

- (o.i2) 


Wt  r ppt [-  (t‘-e)JL,Z*^(  tt)  f H,1  (*'-*) 

- K,  an  ($  -<t>)  j _ 


(0.13) 


V-  Pi 


U,  ~ PR,  [-K,  Q, + R* j£*«  aPk  ($>-$) 


2 ,10'-  Pi 

- K,  **»(*-*)]p^ 


(0.14) 
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These  expressions  for  u-j , u^,  and  are  valid  for  , ip^,  <J>-|  and 

7T 

4>2  less  than  2-  In  order  to  extend  the  domain  of  the  integration  to 
it,  symmetry  of  the  displacement  fields  must  be  considered.  For  , 

TT 

<f>1  and  ^ less  than  ^ these  conditions  may  be  expressed  as 


(0.15) 

U,  (YY/VJ,  P)-  - U,  (T,  TT-  P) 

(0.16) 

Wt<£Mi>P)=  r) 

(0.17) 

(£»  <&A,P)  - -UiMj* P) 

(0.18) 

The  above  analysis  is  valid  only  for  ip  < ip^  < ip 2-  If 
^ > ^2  > ^1  (or  (p  > $2  > <l>i ) then  the  absolute  values  of  ip  - ip. 

(i  = 1,  2)  and  <J>  - cf ^ in  equations  (0.11  )-(0.14)  should  be  considered 
in  the  evaluation  of  w.| , w^,  and  u^.  In  addition  to  the  above, 
the  direction  of  u-j  and  U2  must  be  accounted  for  by  considering  these 
quantities  as  being  + u-j  (or  + u2)  when  ip  <_  < ip  (<J>  < <j>.j  < <J>g ) and 

as  - u.|  (-  u2)  when  ip  > ip2  > ip -|  (4>  > 4>-|  > 4>2)- 


APPENDIX  P 


DERIVATION  OF  A PLANAR  APPROXIMATION  TO  THE 


ELEMENTAL  AREA  ON  A SPHERICAL  SURFACE 


Consider  the  area  ABCD  on  a spherical  surface  of  radius  R as 
shown  in  figure  P.l.  It  is  desired  to  approximate  this  curved  area 
with  a plane  element  so  that  the  integral  of  the  Boussinesq  influence 
functions  can  be  approximated  for  the  spherical  element. 

The  area  on  the  sphere  is  bounded  by  ± A and  6 ' ± A where 
A « 1 radian.  A cone  generator  is  defined  such  that  its  apex  N falls 
on  the  n axis  and  points  A,  B,  C,  and  D fall  on  the  surface  of  the 
cone  (see  fig.  P.l).  The  surface  ABCD  on  the  cone  closely  approxi- 
mates the  surface  on  the  sphere  for  small  A.  Let  point  G be  located 
midway  between  A and  D,  and  let  G'  be  its  projection  on  the  n axis 
(as  shown  in  fig.  P.2).  From  figures  P.l  and  P.2,  the  following  rela- 
tions may  be  established: 


AG  = GJ)  - R juhi  A. 

06  - R<on  & (P- 

NG  = OG  to^L  T‘  9 RsCOQ-A&ulY' 

NA  - NG-AG  = - '**«*)* 

Nt>  - NG  +AG 

_ ^ sooty*' 


G 1 G (% &} 


- ^ xuxl  'Y 


NG 
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BASE  OF  CONE 
GENERATOR 


Fig.  P.l.  Elemental  area  on  spherical  surface 
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Fig.  P.2.  Section  view  of  cone  generator 


Fig.  P.3.  Annular  segment 
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Finally,  for  small  A,  the  area  ABCD  on  the  cone  can  be  closely  approx- 
imated by  a planar  annular  element  including  points  N,  A,  B,  C,  and  D. 
Figure  P.3  illustrates  an  annular  segment  ABCD  with  inner  radius  r-j  , 
outer  radius  r^  and  half  angle  of  p.  In  relating  the  equation  (P.l) 
to  the  geometry  of  the  annular  element  the  quantities  r-j , r^  and 
can  be  described  as 


f\  - NA  = R Aim  ( Y /mhl  Y ' 

(P.2) 

fi  -A/Z)  r R 4Am,  (Y'  + a)  / ao4  nr ' 

(P.3) 

J*  - A 404  Y*  ' 

(P.4) 

The  corresponding  areas  of  the  original  spherical  element  may 
be  compared  to  the  area  of  the  planar  annular  element.  For  the  sphere 


TdY 


'a'-* 


= **2*f-^r#* 

- 4 A 

For  the  annular  element 

A = v'(  r^-r,1) 

/ -R\ 4vnl(Y-A) 


(P.5) 


" 2ir 

- zi  /C#i  Y' 


/tonz  Y' 


(P.6) 


For  small  A,Aa  * As- 


NOMENCLATURE 


A, 


U 


A 


A 


'2  k 


B.  . 
ij 


B 


n 


2k 


U 


do 

E 

E. 

l 

E' 


area  of  annular  segment 
area  of  cell  i 
area  of  cell  ij 

coefficients  defined  in  equations  (4.20)  and  (A. 3) 

area  of  sector  on  sphere 

coefficients  in  equations  (4.18)  and  (A.l) 

radius  of  contact  region  in  chapter  5, 
semi -major  axis  of  ellipse  in  appendix  D 

non-dimensional  radius  of  contact  region 

coefficients  defined  by  equation  (2.12) 

coefficients  defined  by  equations  (4.21)  and  (A. 4) 

coefficients  in  equations  (4.19)  and  (A. 2) 

coefficients  defined  by  equation  (2.10) 

half  width  of  contact  region  in  appendix  K, 
semi -minor  axis  of  ellipse  in  appendix  D 

distance  between  field  point  1 and  a point  in  cell  s.. 

* J 

distance  between  field  point  1 and  centroid  of  cell  s.. 

* J 

hypothetical  interpenetration 

initial  interpenetration  of  spheres 

Young' s modul us 

Young's  modulus  of  body  i 

equivalent  Young's  modulus  for  plane  stress 
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E 

equivalent  Young's  modulus  for  plane  strain 

e 

difference  in  radii  of  curvature,  - R-j 

F 

force  in  context  of  point  loads,  force  per  unit 
length  in  context  of  line  loads 

Fi 

* 

F 

force  at  node  i 
nondimensional  force 

fconsphere 

force  resulting  from  analysis  by  CONSPHERE 

FHertz 

force  resulting  from  Hertzian  analysis 

* 

fconsphere 

nondimensional  force  resulting  from  analysis  by 
CONSPHERE 

F*Hertz 

nondimensional  force  resulting  from  Hertzian  analysis 

f 

profile  function  or  initial  separation  between  two 
points  which  merge  after  deformation 

a 

f 

vector  describing  initial  spearation  fo  two  points 
which  merge  after  deformation 

f. 

i 

initial  separation  of  field  points  i 

f\ 

i 

fN+l  " fi 

r t w 

components  of  f in  r,  t and  w directions  respectively 

f(x,y) 

f evaluated  at  point  (x,y) 

f.j(x,y) 

distance  of  surface  i from  x-y  plane  evaluated  at  (x,y) 

G 

modulus  of  rigidity 

G ( ) 

influence  function  for  displacements  normal  to  a 
surface  in  a cartesian  coordinate  system 

G ( ) 

influence  function  for  displacements  normal  to  a 
surface  in  a polar  or  spherical  coordinate  system 

Gi  ( ) 

G ( ) for  body  i 

g ( ) 

nondimensional  influence  function 
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nondimensional  influence  function  for  displacements  in 
the  radial  direction 

nondimensional  influence  function  for  displacements 
in  the  tangential  direction 

influence  function  for  displacements  tangential  to  a 
surface  in  a cartesian  coordinate  system 

influence  function  for  displacements  tangential  to  a 
surface  in  a polar  or  spherical  coordinate  system 

H(  ) for  body  i 

numerically  generated  nondimensional  influence  func- 
tion for  displacements  in  the  radial  direction 


numerically  generated  nondimensional  influence  func- 
tion for  displacements  in  the  tangential  direction 


total  integral  over  an  annular  segment 


dA.  ./C. 
ij  ijl 


elastic  parameter 


TT  E 


elastic  parameter  (ItH). 


elastic  parameter 


1+v 

ttE 


K.,  evaluated  for  body  j 

1-v2  1-v2 

elastic  parameter  ^ — 

initial  value  of  constant  in  point  mating  procedure 

th 

constant  in  point  mating  procedure  for  j solution 
on  body  i 

inverse  of  Poisson's  ratio  in  chapter  4 and  appendix  A 
constant  in  mapping  function  in  appendix  D, 
number  of  equal  sectors  in  ft  in  chapter  7 


N 

number  of  cells  in  simply  discretized  solution 

Ncon 

A 

n. 

l 

number  of  cells  within  r 

con 

unit  normal  vector  to  surface  of  body  i 

P 

constant  pressure 

pi 

constant  pressure  in  cell  i 

P 

max 

* 

P. 

1 

★ 

P 

max 

maximum  pressure 

nondimensional  pressure,  kP^ 
★ 

kP 

max 

P2k  (cos  e) 

Legendre  polynomial  in  cos  0 

P ( ) 

interfacial  contact  pressure  in  polar  or  spherical 
coordinate  system 

P(x,y,z) 

interfacial  contact  pressure  in  cartesian  coordinates 

R 

radius  of  curvature 

Ri 

radius  of  curvature  of  body  i (k  = 1,2) 

2R,  R« 

Rm 

mean  radius  of  curvature,  D-  . D- 

R2  K1 

(note:  in  conformal  theory  -R0  is  substituted  for 

r2)  2 

AR 

R2  " R1 

r 

coordinate  of  generic  point 

* 

r 

nondimensional  coordinate  of  generic  point,  r/Rm 

A 

r 

unit  vector  in  mean  radial  direction 

r 

position  vector 

.» 

r. 

l 

position  vector  to  point  on  surface  of  body  i 

rb 

radius  of  blending  point  of  pit  on  sphere 
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r 


c 


r 


i 


rA<(> 


radius  of  curvature  at  "edge"  of  pit  in  chapter  7, 
radius  of  centroid  in  appendix  G 

radius  used  in  convergence  studies 

radius  of  inner  boundary  of  contact  region 

chords  used  in  analysis  of  influence  functions  for  a 
cylinder  in  chapter  4 and  appendix  B (i  = 1,2), 
radius  of  inner  boundary  of  ith  cell  in  chapter  7 


radius  of  outer  boundary  of  contact  region 

x and  z components  respectively  of  vector  r. 

r/Rm 

m 

r./R 
b m 

r /R 
c m 

length  of  cell 
width  of  cell 


SCF 


TOL 

Tl>  V T3 


t 


separation  of  merging  field  points 
separation  vector  of  merging  field  points 
r,  t,  and  ^ components  of  § 

singularities  in  influence  function  for  radial  dis- 
placements on  a sphere 

stress  concentration  factor 

region  included  in  ring  i between  rays  j and  j + 1 
distance  along  contour  curve  of  body  i (i  = 1,2) 
si  on  jth  solution 

tolerance  used  to  determine  merging  of  field  points 

singularities  in  influence  function  for  tangential 
displacements  on  a sphere 

unit  vector  in  mean  tangential  directions 
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u 

u 


1 

vj 


V . 
1 


w 


w. 

1 


w. 


w. 

1 


(x,y,z) 


Z 


Z 


0 


Z 


1 


Z 


1 


z . 

l 


a 


ot-|  ,a2 
3 


displacement  tangent  to  surface 
displacement  in  radial  direction 
displacement  in  tangential  direction 
displacement  in  z direction 
displacement  in  tangential  direction 
tangential  displacement  on  body  i (i  = 1,2) 
coefficients  defined  in  equation  (2.11) 
displacement  in  w direction  on  body  i (i  = 1,2) 
displacement  normal  to  surface 
normal  displacement  on  body  i (i  = 1,2) 
w.j  at  outer  boundary  of  contact  region 
nondimensional  normal  displacement  w./R. 

K l i 

coordinates  of  generic  point  in  cartesian  coordinate 
system 

unit  vectors  on  body  i 

i 0 

generic  vector  in  Z plane.  Re 

A 

vector  in  x direction  in  Z plane 
vector  in  Z plane 
conjugate  of  Z-j 
val ue  of  f . (x,y) 

angle  between  mean  radial  direction  and  the  Z axis 

angles  defined  in  figures  4.12  and  B.l 

coordinate  of  generic  point  in  spherical  coordinate 
system 

centroidal  radius  of  sector  S. . 

* J 
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Y 

▲ 

angle  between  F and  Z axis 

A A 

Vy2 

unit  vectors  in  cartesian  coordinate  systems  fixed  to 
bodies  1 and  2 respectively 

Yxy 

shear  stress 

A' 

small  angle  «1  radian 

A 

constant  defined  in  equations  (4.22)  and  (A. 5) 

Ai 

rigid  body  translation  of  body  i (i  = 1,2), 
triangle  i (i  = 1,6)  defined  in  appendix  G only 

6 

approach 

★ 

6 

6/  R 

m 

6consphere 

* 

6 CONSPHERE 

approach  predicted  by  CONSPHERE 
6CONSPHERE/Rm 

6Hertz 

* 

6 Hertz 

approach  predicted  by  Hertzian  theory 
^Hertz^m 

6 9 6. 

r t 

approach  in  the  r and  t directions  respectively 

e 

root  mean  square  error 

£ , £ 

x y 

components  of  strain  in  appendix  I 

£ . . 

1J 

strain  tensor  (i,j  = 6,  ip,  r,  u>,  $) 

c 

generic  vector  in  c-plane 

A A 

5r  c2 

unit  vectors  in  cartesian  coordinate  system  fixed  on 
bodies  1 and  2 respectively 

A /\ 

nrn2 

unit  vectors  in  cartesian  coordinate  system  fixed  on 
bodies  1 and  2 respectively 

0 

generic  angle  used  in  spherical  and  polar  coordinate 
systems 

K 

elastic  parameter  equal  to  ^ for  plane  stress  and 
3-4v  for  plane  strain 
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angle  between  normal  vector  n.  and  r direction 
/ 1 

modulus  of  rigidity,  Lame's  constant 
Poisson's  ratio 

Poisson's  ratio  of  body  i (i  = 1,2) 
equivalent  Poisson's  ratio  for  plane  stress 
equivalent  Poisson's  ratio  for  plane  strain 

-1  /V 

angle  between  f and  r 

half  vertex  angle  of  annular  segment  in  chapter  5 
and  appendix  P, 

a boundary  radius  (rT  or  r ) of  the  contact  region  in 
chapter  7 1 0 

components  of  stress  in  appendix  I 

vector  in  c-plane 

conjugate  of  o-j 

angle  defined  in  equation  (7.9)  in  chapter  7, 
angle  describing  orientation  of  u^  in  appendix  N 

shear  stress 

regularization  parameter 

generic  angle  in  spherical  and  polar  coordinate  sys- 
tems fixed  to  body  2 

angle  defining  boundary  of  contact  region  on  body  2 
27T/m  in  chapter  7 

generic  angle  in  spherical  and  polar  coordinate 
systems  fixed  to  body  1 

angle  defining  boundary  of  contact  region  on  body  1 

function  which  is  minimized  in  Functional 
Regularization  method 


contact  region 


projection  of  contact  region  onto  x-y  plane  in 
chapters  2 and  7 only 

region  of  cell  i in  chapters  2 and  7, 
contact  region  on  body  i in  chapters  3-6 

contact  region  of  cell  j on  body  i 

candidate  or  tentative  contact  region 

tentative  contact  region  of  cell  j on  body  2 

unit  vector  defined  by  equation  (3.3b) 

mapping  function 

j A 

angle  between  r^  and  x (i  = 1,2) 
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